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ABSTRACT 

A  two-fluid  collislonless  model  of  a  plasma 
is  used  to  describe  the  behavior  of  a  fully  ionized  gas 
under  axially  symmetric  radial  compression^  a  problem 
referred  to  as  the  pinch.   Two  pressure  tensors  are  derived 
and  used  in  the  equations,  one  for  the  ions  and  one  for  the 
electrons.   The  ion  pressure  is  determined  by  taking  higher 
moments  of  the  ion  Vlasov  equation  and  assuming  that  the 
third  moments  of  the  distribution  function  can  be  neglected. 
Expressions  for  the  electron  pressure  are  found  from  an 
asymptotic  expansion  of  the  electron  distribution  function. 
The  expansion  parameter  is  the  square  root  of  the  ratio  of 
the  electron  mass  to  the  ion  mass.   The  electron  pressure 
tensor  is  found  to  depend  on  the  magnetic  fields  as  well  as 
the  density  while  the  ion  tensor  is  approximately  propor- 
tional to  the  density  raised  to  a  power  which  depends  on 
the  radius. 

The  equations  of  motion  are  put  into  Lagrangian 
coordinates  and  solved  numerically  on  the  GDC  6600  digital 
computer.   A  proof  of  the  stability  of  the  difference 
method  in  the  special  case  of  zero  electron  pressure  is 
included . 
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A   few  types  of  pinches  are  considered.   The  first 
is  a  theta  pinch^  deriving  its  name  from  the  fact  that  the 
current  used  to  compress  the  plasma  is  in  the  circumferen- 
tial direction.   In  this  pinch^  and  in  the  others  to  follow, 
we  observe  the  familiar  magnetic  compression  waves  associated 
with  the  two-fluid  model.   A  train  of  waves  is  generated  at 
the  boundary  and  propagates  inward  as  the  radius  of  the 
plasma  column  decreases.   The  pulses  peak  strongly  upon 
reaching  the  center  and  are  then  reflected  outward  whereupon 
they  interfere  with  and  pass  through  the  remaining  incoming 
tail  of  the  wave  train.   In  the  theta  pinch  the  radius  of 
the  column  reaches  a  minimum  value  after  the  applied 
magnetic  field  attains  its  final  value.   It  then  oscillates 
about  a  radius  somewhat  below  its  initial  position. 

The  second  pinch  considered  is  a  hybrid  z-theta  pinch 
in  which  a  plasma  which  contains  a  constant  axial  magnetic 
field  is  subjected  to  a  rising  circumferential  field  on  the 
boundary.   As  the  plasma  column  contracts  an  increasing 
axial  field  is  applied  at  the  boundary  in  such  a  manner  as 
to  preserve  the  initial  axial  magnetic  flux  in  the  gas.   It 
is  found  that  in  this  way  very  strong  compressions  can  be 
obtained  without  the  solutions  becoming  discontinuous. 

The  last  flows  considered  are  produced  by  compressing 
a  plasma  which  initially  contains  both  components  of  the 
magnetic  field  in  a  force-free  configuration.   To  initiate 
the  collapse  one  component  of  the  field  is  increased  on  the 
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boundary  while  the  other  Is  held  fixed  there.   This 
results  in  some  unusual  local  behavior.   Depending  on  the 
initial  and  boundary  conditions  the  density  either  rises 
sharply  or  falls  strongly  at  the  boundary  as  the  plasma 
is  compressed. 

In  all  cases  the  compressions  are  sufficiently  weak 
so  that  the  solutions  are  well  behaved^  that  is^  no  shocks 
develop.   Consequently ,  there  is  considerable  penetration 
of  the  magnetic  fields  at  the  boundary  into  the  plasma. 
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CHAPTER  I.   INTRODUCTION  AND  SUMMARY 

The  study  of  the  "behavior  of  a  plasma  under 
cylindrical  compression,  the  pinch  problem,  received  its 
original  stimulus  from  studies  of  possible  methods  to 
control  and  contain  thermo-nuclear  fusion.   The  problem  is 
still  under  active  investigation,  both  theoretically  and 
experimentally.   The  term  pinch  is  used  to  describe  the 
motion  of  a  plasma  in  either  a  cylindrical  or  toroidal 
configuration  when  an  increasing  magnetic  field  is  applied 
to  its  boundary.   Since,  for  an  ionized  gas,  a  magnetic 
field  behaves  like  a  pressure  the  gas  is  accelerated  inward 
by  the  field.   In  practice  the  fields  are  produced  by  either 
running  a  current  through  the  plasma  or  through  wires 
surrounding  the  plasma.   The  type  of  current  used  gives 
the  pinch  its  name;  when  the  current  is  in  the  circumferen- 
tial direction  we  speak  of  a  theta  pinch  while  an  axial 
current  gives  rise  to  a  z  pinch. 

The  models  used  in  the  theoretical  description  of 
the  pinch  vary  in  complexity  from  the  one-fluid  infinite 
conductivity  magnetohydrodynamic  model  represented  by  the 
Lundquist  equations  to  self -consistent  particle  theories 
involving  Vlasov  equations  for  each  component  of  the  gas. 
Our  treatment  uses  a  two-fluid  collision-free  model  which 
is  almost  consistent  with  the  particle  formulation. 

Much  of  the  research  on  the  pinch  problem  has 
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involved  numerical  studies  of  the  motion.   One  of  the 
earliest  was  performed  by  Rosenhluth  (19)  using  his  so- 
called  snow  plow  model  in  which  all  material  swept  up  by 
the  incoming  boundary  is  piled  in  a  thin  layer  and  travels 
with  the  boundary.   In  the  same  paper  he  also  investigates 
numerically  a  simplified  particle  model.   This  model  predicts 
the  formation  of  shocks  which  precede  the  incoming  boundary 
and  are  reflected  from  the  center. 

For  a  more  accurate  description  of  an  ionized  gas 
a  two-fluid  model  such  as  was  proposed  by  Gardner  et  al.  (7) 
is  used,  that  is  the  ions  and  electrons  are  assumed  to  be 
present  in  two  distinct  gases.   In  a  pressure-free,  steady 
state,  one  dimensional  solution  using  this  model  Adlam  and 
Allen  (1)  found  flows  containing  periodic  oscillations  or 
solitary  waves.   Their  work  was  extended  by  Banos  and  Vernon 
(5)  who  included  an  isotropic  fluid  pressure  and  by  Hain, 
Ldst,  and  Schliliter  (ll)  who  considered  also  non-isotropic 
pressures  derived  from  the  Vlasov  equations.   Gardiner 
and  Morikawa  (&)  examined  the  effect  of  temperature  on  these 
solitary  waves  by  expanding  the  pressure  in  an  asymptotic 
series  in  a  small  parameter  representing  the  difference 
between  the  fluid  velocity  and  the  Alfven  Mach  speed.   Auer, 
Hurwitz,  and  Kilb  (2),  Rossow  (20),  and  Morton  (I5,l6),  have 
performed  numerical  simulations  of  the  magnetic  piston  problem 
using  two-fluid  models.   Morton  related  the  time  dependent 


solution  to  the  steady  state  results.   He  took  the  pressure 
to  "be  that  of  an  isentropic  fluid  with  y  equal  two.   Although 
this  does  not  give  an  accurate  description  of  a  warm  plasma 
it  is  included  to  allow  a  qualitative  description  of  the 
flow  after  breaking  has  occurred.   In  the  case  of  zero 
pressure  his  model  is  consistent  with  a  particle  formulation. 
Our  pressure^  however^  is  derived  from  the  Vlasov  equations 
and  is  therefore  more  nearly  consistent  with  the  particle 
formulation  for  warm  plasmas. 

Numerical  simulations  of  the  pinch  have  also  been 
performed  using  two-fluid  models.   Kilb  (12) ;,  using  such  a 
representation,  studied  the  theta  pinch  assuming  that  the 
electron  pressure  was  proportional  to  the  square  of  the 
magnetic  field  while  ignoring  the  ion  pressure.   Hain,  Hain, 
Roberts,  Roberts,  and  Koppendorfer  (lO)  included  electrical 
and  thermal  conductivity  in  their  model.   They  assumed  that 
the  electrons  are  heated  ohmically  while  the  ions  are  heated 
by  shocks.   The  pressures  of  the  two  components  are  taken 
to  be  equal  to  the  products  of  their  respective  densities 
and  temperatures.   To  allow  numerical  computation  through 
shocks  and  to  determine  the  ion  heating  an  artificial  vis- 
cosity as  discussed  by  Richtmyer  (l8)  is  included.   Rather 
than  relate  the  current  to  the  fluid  velocities  they  use  a 
tensor  form  of  Ohm's  law.   Following  Morton,  we  express  the 
current  as  a  weighted  sum  of  velocities  of  the  components 
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so  no  empirical  form  of  Ohm's  law  is  needed.   We  assume 
the  plasma  is  collision-free,  so  there  is  neither  heat 
conduction  nor  resistivity. 

Aside  from  the  obvious  effects  of  cylindrical 
convergence,  the  "basic  difference  between  the  one  dimensional 
piston  problem  and  the  axially  symmetric  two  dimensional 
pinch  problem  is  that  since  div  B,  which  in  cylindrical 
coordinates  is  d(rB-,)/dr,  must  vanish  rB-,  is  a  function 
only  of  time.   For  B-,  to  be  regular  in  the  region  considered 
it  must  vanish  identically.   Since  the  fields  are  then 
transverse  to  the  flow  the  scale  length  which  should  be 
used  is  given  by  the  mean  gyro-magnetic  radius  of  the  ions 
and  electrons,  for  this  is  the  scale  length  of  the  transverse 
waves  obtained  by  Morton.   The  geometry  also  differentiates 
between  the  two  transverse  fields;  the  influence  of  the 
axial  field  is  quite  different  from  that  of  the  circumfer- 
ential field.   There  are  two  main  reasons  for  this.   First, 
the  circumferential  field  is  required  to  vanish  at  the 
origin  while  the  axial  field  is  required  only  to  be  regular, 
and  secondly  the  centrifugal  effects  are  expressed  in  terms 
of  the  axial  field  in  the  equations  of  motion. 

In  this  work  an  infinite  cylinder  of  plasma  composed 
of  two  oppositely  charged  fluids,  the  ions  and  electrons,  is 
taken  to  be  initially  at  equilibrium  with  constant  density 
and  constant  axial  current.   This  determines  the  pressure 
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and  magnetic  fields.   For  the  special  case  in  which  the 

pressure  is  zero  it  is  necessary  to  include  a  circumferential 

current  to  maintain  the  equilibrium.   Pressure  effects  at 

the  boundary  are  neglected  so  the  edge  appears  as  a  sharp 

division  between  plasma  and  vacuum.   The  flows  to  be  studied 

are  produced  by  increasing  the  magnetic  field  at  the  boundary, 

compressing  the  plasma.   The  nature  of  the  flow  depends  on 

R  ,   the  initial  radius  of  the  column, 

a   and  a  ,   the  initial  values  of  the  ion  and 

electron  pressures  as  functions  of  position,  and 

B^(T,Y    )  and  B^(T,Y    ),   the  magnetic  fields 
2^  ■'  max'      ;>    max'^^        ^ 

imposed  at  the  boundary  as  functions  of  time.   Also  occurring 

is  the  small  parameter 

2      e'^m"  ±      t 

e,   defined  by  e  ^  -  —  ,  where  e   and  rri  are 

e  m 

the  charges  and  masses  of  the  particles  in  the  two  components. 

In  Chapter  II  the  equations  of  motion  are  derived 

from  collision-free  Vlasov  equations  for  each  component. 

The  system  is  completed  by  using  Maxwell's  equations  without 

displacement  current,  thus  restricting  ourselves  to  studying 

phenomena  which  occur  on  a  length  scale  which  is  large 

compared  to  the  Debye  length.   Cylindrical  coordinates  are 

used  and  the  solution  is  assumed  to  depend  only  on  time  and 

radial  distance.   In  standard  fashion  the  first  two  moment 

equations  yield  conservation  equations  which  contain  mean 

velocities.   They  also  contain  pressure  terms.   To  determine 

the  ion  pressure  we  take  higher  moments  of  the  ion  Vlasov 
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equation  and  assume  some  symnietry  in  the  radial  velocity 
distribution.   This  allows  us  to  prove  that  the  third  moment 
terms  are  negligible.   Alternatively ^  we  could  have  assumed 
low  ion  temperatures  to  arrive  at  the  same  result.   The  ion 
pressure  is  found  to  "be  strongly  influenced  by  the  cylindri- 
cal geometry^  with  the  diagonal  components  of  the  pressure 
tensor  being  nearly  proportional  to  the  density  raised  to 
a  power  which  depends  on  the  radius.   The  electron  pressure 
is  also  proportional  to  the  density,  but,  unlike  the  ion 
pressure,  the  factor  of  proportionality  depends  on  the 
magnetic  fields.   The  components  of  this  tensor  are  found 
by  an  asymptotic  expansion  in  e  of  the  electron  distribution 
function  in  a  manner  suggested  by  Morawetz  (l4). 

The  ion  and  electron  equations  are  put  into  Langrang- 
ian  coordinates,  retaining  only  the  lowest  order  terms  in  £. 
Taking  initial  conditions  into  account  the  system  is  reduced 
to  three  equations  for  R,  B^  and  B^  where  R  is  the  radius 
expressed  as  a  function  of  the  new  coordinates  and  B  and 
B^  are  the  circumferential  and  axial  magnetic  fields.   The 
equations  are  of  mixed  hyperbolic-parabolic  type  with  the 
only  non-trivial  characteristics  being  the  forward  and 
backward  sound  speeds. 

The  one-fluid  analog  to  the  problem  is  briefly 
examined.   It  is  shown  that,  with  proper  scaling  of  the 
independent  variables,  the  one-fluid  solution  is  obtained 
from  the  two-fluid  one  by  a  singular  limiting  process  in 
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which  the  scaling  parameter  goes  to  zero. 

Chapter  III  contains  the  numerical  methods  used  to 
solve  the  equations  derived  in  the  preceding  chapter.  A 
proof  of  the  stability  of  the  method  used  is  given  for  the 
special  case  in  which  the  electron  pressure  is  absent  and 
the  ion  pressure  is  initially  constant.   The  proof  depends 
on  the  energy  method  and  uses  an  energy  integral  of  the 
system. 

The  last  chapter  contains  the  results  of  the  numeri- 
cal simulation  of  the  two  types  of  pinch.   For  the  theta 
pinch  we  observe  the  characteristic  two-fluid  compression 
waves  being  generated  at  the  outer  boundary  and  propagating 
inward  with  a  speed  somewhat  in  excess  of  the  Alfven-sound 
speed  predicted  by  the  linearized  theory.   The  speed  depends 
on  the  amplitude  of  the  waves  which  can  become  quite  large 
without  initiating  the  formation  of  shocks.   When  the 
disturbance  reaches  the  center  the  density  rises  very 
strongly.   As  each  wave  arrives  at  the  center  it  is  reflected 
and  travels  back  through  the  plasma  interacting  with  the 
remaining  tail  of  the  wave  train. 

We  investigate  also  a  hybrid  z-theta  pinch  in  which 
a  plasma  containing  a  constant  axial  magnetic  field  at 
equilibrium  is  compressed  by  imposing  an  increasing  circum- 
ferential field  at  the  boundary.   We  impose  also  the  condition 

that  B^   vanishes  at  the  boundary.   This  is  almost  equivalent 
^r 

to  the  requirement  that  the  axial  magnetic  flux  in  the  plasma 
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is  conserved.   With  these  conditions  we  find  that  very  strong 
compressions  with  large  density  rises  are  produced.   What  is 
more,  the  solution  remains  well  behaved  for  a  time  far  longer 
than  one  would  expect  from  the  results  of  the  theta  pinch 
calculations . 

Finally,  we  look  at  some  flows  which,  while  they  might 
not  be  experimentally  realistic,  exhibit  some  rather  unusual 
phenomena .   In  these  flows  we  have  both  magnetic  fields 
present  initially  in  a  force  free  configuration.   By  specify- 
ing that  the  initial  axial  current  is  constant  the  state  is 
uniquely  determined.   We  then  raise  one  of  the  fields  at  the 
boundary  and  require  that  the  other  field  is  held  fixed  there 
at  its  initial  value.   By  imposing  this  condition  we  observe 
some  paradoxical  effects.   When  Bp  is  raised  the  plasma 
develops  a  sharp  rise  in  density  at  the  boundary  while,  for 
a  particular  initial  configuration,  raising  B^  produces  the 
opposite  effect,  the  density  at  the  boundary  falls.   These 
results  are  explained  by  the  effects  of  the  initial  field 
geometry  and  the  special  nature  of  the  boundary  conditions. 
The  compression  waves  produced  in  these  cases  are  similar  to, 
but  weaker  than,  those  observed  in  the  theta  pinch. 
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CHAPTER  II.   THE  EQUATIONS  OF  MOTION 

In  this  chapter  we  discuss  the  model  used  to  represent 
the  plasma  and  derive  expressions  for  the  pressure  tensors 
which  appear  in  the  equations  of  motion.   Cylindrical 
Lagrangian  coordinates  are  used  to  put  the  equations  into 
a  form  suitable  for  numerical  computation. 

a .   The  two-fluid  collisionless  model 

We  consider  a  model  in  which  the  plasma  is  composed 
of  two  oppositely  charged  components*  the  dynamics  of  each 
being  governed  by  collisionless  Vlasov  equations  for  their 
distribution  functions  f~(tjr,u).   The  plus  and  minus 
superscripts  denote  ion  and  electron  terms.   The  arguments 
of  these  distribution  functions  are  time^  the  position 
vector^  and  the  velocity  vector.   The  electric  and  magnetic 
fields,  E(tjr)  and  B(t,r),  and  the  current  j(tjr)  obey 
pre-Maxwell  equations.   The  neglect  of  the  displacement 
current  restricts  us  to  studying  low  frequency  phenomena. 
We  have  then,  in  rationalized  mks  units. 
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The  new  dimensionless  variables  will  be  denoted  by  the  same 
symbols . 

Since  the  pinch  problem  deals  with  the  radial 
compression  of  the  plasma  it  is  natural  to  assume  that  the 
solution  is  axially  symmetric,  being  dependent  only  on  the 
radius  "r"  and  time  "t".   Let  the  subscripts  1,  2  and  3 
denote  the  r,  9   and  z  components  of  a  vector.   The  force 
terms  can  then  be  written 

F^  =  Ke^lE^  +  u^B^-  u^B^)  +  |  u^ 

(2.5)      ^2  =  Ke^(E2  +u^B^  -  u^B^)  -  i-  u^u^ 

F^  =  Ke^(E^  +u^B2-  u^B^) 

where 

1   for  the  +  equation        ^      +  - 
K  =  <^  and  e  =  -      . 

-1   for  the  -  equation  e  m 

The  force  components  are  seen  to  be  composed  of  E  +  u  x  b 

terms  plus  centrifugal  effects  arising  from  the  cylindrical 

2 
geometry.   For  singly  ionized  particles  e   is  the  mass  ratio. 

The  first  equation  of  (2.2)  is  (rB-,  )   =  0.   To  avoid 

the  existence  of  a  magnetic  monopole  at  the  origin  B-,  is 

taken  to  be  identically  zero.   Thus,  using  the  force  terms 

from  (2.5)  the  Vlasov  equations  take  the  form 

2 

fl   +  U,  f^   +  -^  rU  ^   ^   ----      .   r.„^r/^     ...   ^    ..   ^   N^ 


(2.6) 


t  '"  "l"r  ^  r"  "  ^1 r~  ^^2  "^  ^^  ^^\   +U2B^-u^B2)f^u^ 


+  (Eg-  U-|_B  )f^U2  +  (E  +u^B2)f^u  ]  =  0 
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The  macroscopic  equations  of  motion  are  obtained 

from  (2.6)  by  taking  moments.   We  define  first  the  ion 

and  electron  number  densities  and  their  mean  velocities! 

/  r 

rr  =      f^   ldu|     ,     [XT  ,\r  ,w )  =  \  \x  r   Iduj  /  M 


Direct  integration  of  (2.6)  gives  the  conservation  of 
mass  equations: 

(2.7)  (^N^)t  +  [r^^)^   -  0 

The  assumption  of  charge  neutrality  gives  N  =  N~,  so 

subtracting  the  minus  from  the  plus  equation  in  (2.7)  yields 

(r(U"'"  -  U~) )   =  0.   Again^  to  avoid  singularities  at  the 

origin,  vje  must  have  U  =  U~ .   We  denote  these  common  values 

by  N  and  U.   Momentum  conservation  equations  are  obtained  by 

multiplying  equations  (2.6)  by  u-,,  Up  and  u  and  integrating, 

assuming  that  the  distribution  functions  vanish  sufficiently 

rapidly  as  the  velocity  goes  to  infinity.   This  gives 

2      2 
(rNU)^  +  (rNU^)^  -  WT       +  ^^^   i^^Al^ r    '    ^2  ^ 

-  Ke%r[E^  +  \rB  -W^Bp]  =  0 

2 
(rNV-)^  +  (rNUV-)^,  +  NUV'  +  if^  t(^Pi2)r  +  ^iV 
(2.8) 

-  Ke%r[E2  -  UB  ]  =  0 

2 
(rNW^)^  +  (rNUW^)^  +  ^f^  ^""^l^) r 

-  Ke^Nr[E  +UB2]  =  0 

Where  P"." .  =  --^r        [u.-  '  u.r"  du  Mru.-   \x  .r    du  ir"  du 
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Equations  (2,2)  take  the  form 

Bp   =  E    ,      rB   =  -  (rE  ) 
"^t  ^r         -'t        ^  ^ 

(2.9)      (rBp)„  =  rJ   =  -^  rN(W+-  W" ) 
B   -  -  Jp  =  -"-%  N(V+- V") 

The  system  (2.7),  (2.8)  and  (2.9)  contains  eleven  equations 

for  the  nineteen  quantities  N,  U,  V-,  W- ,    B^,  B^,    E-^,    E^,    E^, 

+  +  +  ±  ,  .       . 

P7n  ,    Ppp:>  FTp^  P-ix-   1'°  close  the  system  xt  is  necessary  to 

obtain  expressions  for  the  eight  pressure  terms. 

It  is  useful  at  this  point  to  introduce  Lagrangian   '- 

coordinates  which  simplify  the  equations  and  have  the 

desirable  effect  of  leaving  the  boundary  of  the  plasma 

fixed  in  the  new  coordinates.   Let  Y  be  defined  for  each 

point  in  the  plasma  as   l/(2Tr)  times  the  total  number  of  ions 

lying  between  the  origin  and  the  point  and  let  T  be  the  total 

time  derivative.   Thus,  if  we  let  R(TjY)  =  r,  we  have 


R(T,Y) 
Y  = 


ri  Y 

rN(t,r)  dr  ,     5-^=0- 


0 
The  transformation  to  Lagrangian  coordinates  can  be  expressed 

in  the  more  useful  differential  form  as 

a_  _  ^   d         A  -  ^_  4.  TT  ^ 
Sy  "  M  5F  ^       ^T  '  ^^  ^     ' 

from  which  we  readily  see  that 

N  =  ^  ,  U  =  R^   . 
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Equation    (2.7)   becomes   an  identity   in  the  new   coordinates 
and    (2.8)    and    (2.9)    take   the   form 


R, 


V±^    l+s2 


K, 


TT 


R-  +  ^^    [(RPii)y    -    Ry^22]    -    K£^[Ei+V"B^-¥-B2]  =  0 


1 
(2.10)<(RV-)^    +^rK(R^Pl2V    "    Ke^R[E2-  R^B^]    =    0 


Kr 


^T   +  ^^    ^^   ^hh   -   Ke^EE^+R^Bg]    =   0 


r 


"y'^    P  T 


-       RrpB^  y 


E 


3Y 


(2.11)< 


RRyB^rjn    -    RRpB^Y  -    -     v^>-^2''y 


(RBp)v  = 


3Y 


'ojy   ~  2 


RB^Y 


^£_  rv+ 


-  (RE^), 
w"^-  w") 


1+e 


This  system  allows  considerable  simplification. 

The  electric  fields  in  (2.10)  are  eliminated  by  multiplying 

2 
the  electron  equations  by  e  and  adding  them  to  the  corres- 
ponding ion  equations.   This  gives 


2. 


(1+£")R^^    -   I   (V+^+    eV^)    +    (1    +£^)[(R(P^+P--^))y 


-,+ 


■+ 


T+    T.r-' 


(2.12) 


-    Ry(PJ2+P22)]    -    ^B^(V   -V")    +    eB2(W'*--W-)    =    0 
(R(V++   eV))^  +    (l+e^)(R^(P+   -    P"    ))v  =   0 


T 


12    12''''Y 


^+ 


(W++  e^W")^  +  (l+e2)(R(Pj^+  ^l^))^   -  0 


To  eliminate  the  electric  fields  from  the  remaining  equations 
we  subtract  the  electron  equations  from  their  ion  counterparts 
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in  the  second  and  third  memhers  of  (2.10).   The  results  are 
differentiated  with  respect  to  Y  to  yield 

2 

(R(v+-  v-))^Y  +  (^^(^^2-  4  ^I2))yy  -  H^  ((RE2)y-  (RV^)^ 

2 
(W+-  W-)^Y  +  (^(^13  -  ^  ^IJ^^YY  -  ^^  ^^3Y  +  V2)y)  =  ^ 

Using  the  equations  in  (2.11)  we  obtain  the  following 
equations  for  the  magnetic  fields! 

[(RB2)Yy  -  Ry^2]t  +  ~^~2  ^^^^13  "  '2    ^IJ^  ^YY  =  ° 
(2.13) 

[(R^B-^y)y  "  I^V^^T  ■  ~"2  [R^(Pl2  -  ■%  ^12)]yY  =  ° 

The  first  equation  in  (2.12)  can  be  simplified  using  (2.11) 
and  the  identity 

V+^  +  eV^  =  -i-^-  [(V+  +  sV)^  +  e2(V+  -  V-)^]. 
l+e 

It   then  becomes 

o 
Rrpiji    -    R(B^y)       +  B2(RB2)y   +   RB^B^Y 

(s.i'O 

+   (R(P+,   +   P-,))^-    K^(4,   +   P;^)    =  Jjf^j-P- 

¥e  can  proceed  no  further  without  finding  explicit 
representations  for  the  pressure  terms.   It  should  be  noted 
here  that  if  we  assume  that  all  pressure  terms  are  zero, 
the  second  equation  in  (2.12)  shows  that  the  right  hand  side 
of  (2.1-^!)  can  be  taken  to  be  identically  zero.   We  have  then 
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in  (2.13)  and  (2.1-^!)  three  equations  for  R,  B^  and  B^. 

Since  both  the  ion  and  electron  pressure  terms 
will  be  determined  using  series  developments  in  the  small 
parameter  £^  it  is  necessary  to  know  the  order  of  all 
quantities  appearing  in  the  equations.   R,  Bp  and  B^  are 
of  order  one  since  they  were  used  to  define  the  scaling 
parameters.   From  (2.11)  we  find  that  Ep  and  E^  are  also 
of  this  order.   The  first  equation  in  (2.10)  in  which  E-, 
appears  with  R„„  leads  us  to  conjecture  that  E-,  is  0(l/e). 
We  will  use  this  in  determining  the  orders  of  the  remaining 
terms  and  we  shall  see  that  this  conjecture  is  consistent 
with  the  solution.   Define  E  =  e  E-,  .   E  is  0(l)  and  is 
assumed  to  be  regular. 

b.   The  ion  pressure  tensor 

To  obtain  equations  for  the  components  of  the  ion 
pressure  tensor  we  make  the  assumption  that  the  ion  distri- 
bution function  is  symmetric  about  the  mean  radial  velocity, 
that  iSj  f  (t,  r,U^-u-,  ,Up,u  )  =  f  ( t,  r^U-u-,  ,Up,u  ) .   Grad  (8)  and 
Chew,  Goldberger,  and  Low  {U)    have  made  similar  assumptions 
in  cases  where  the  rapid  spiraling  of  the  particles  about 
a  guiding  center  implies  isotropy  of  the  stress  tensor  in 
the  plane  perpendicular  to  the  magnetic  field.   This 
argument  fails,  however,  for  the  relatively  heavy  ions. 
Instead,  we  assume  symmetry  in  the  radial  direction  and 
show  that  this  assumption  implies  symmetry  to  first  order 
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in  e  in  the  other  directions. 

To  prove  that  the  distribution  function  is  symmetric 
in  the  9   and  z  components  of  velocity  to  first  order  in  e 
we  define  two  functions  G  and  H  by 

G{t,T,u^,u^,u^)    =  f"*'(t,r,u^,U2,u^)  -  f'^(t,r,u^,-U2,u^)   and 

H(t,r,u^,U2,u^)  =  f'^(t,r,u^,U2,u^)  -  f''"(t,r,u^,U2,-u^)  . 

G  and  H  measure  the  symiietry  of  the  ion  distribution  function 

in  the  planes  perpendicular  to  the  Up  and  u^  axes.   From  (2.6) 

we  find  differential  equations  satisfied  by  these  functions! 

Up  u  u 

^t  +  ^l^r  +  -  %  +  %  -  — --  ^u^  =  O(^) 

H,  +  u^H^  +  -^  H^^  +  EH^_^  -  -p-   H^^  =  0(s)  . 

If  we  take  G  and  H  to  be  initially  zero  then  it  is  clear  that 
they  will  remain  0(e)  for  all  finite  time.   Thus,  the  ion 
distribution  function  is  symmetric  to  0(e)  about  the  values 
Ur.  -   0   and  u^  =  0.   We  have  then  the  desired  result  that  all 
third  moments  vanish  to  this  order. 

There  is  another  justification  for  ignoring  third 
moments  which  involves  a  temperature  related  parameter. 
We  define  uj(t,r),  vj(t,r)  and  wj(t,r)  by 

Uj   +  Ujuj  +  -^  +  E  =  0  ,   (rvj)^  +  Uj(rvj)^  =  0  , 
t       r 

wj  +  ujwj   =  0  , 

t       r 
and  Pj  q  and  s  by 

-17- 


u^  -  uj  =  ap  ,    u^  -  vj  =  aq  ,    u^  -  wj  =  as 

where  a  is  a  parameter  related  to  the  temperature.   We 

write  (2.6)  in  the  form 

2 
fj  +  u-^f^  +  -|  f+u^  -  -"-^   f+u^  +  E  f+u^  =  0(e) 

and  substitute  into  this  the  new  variables  p^  q  and  s  giving 

a-l[(U+  +  Ujuj  +  -I-  +  E)?p  -  i  (rvj  +  Uj(rvj)^)?^l 


+ 

^0 

t  "^  ""O^r  "^  "   r   "p    r  ^^"0  "  ^"O^^q 


+  f,  +  U+f^  +  2  — -^  f   -  i  (qU+  +  pv;^)f 


+  a[pf^  +^f^-^f^    ^    0(e)  , 


where  f(t,r,p^q^s)  =  f  ( t  ^  r^u-,  ^u^^u  )  .   If  we  now  consider  a 
a  small  parameter  and  express  f  as  a  power  series  in  itj  we 
would  find  that  the  zero  order  term  satisfies 

fo   +  Ujf^  +  2  -^  fp  -  -  (quj  +  pvj)fq  =  0(e)  . 

Using  this  we  find  that  all  third  moments  are  0(  a-^)  and 
can  be  neglected.   The  condition  that  a  is  small  is  that 
the  temperature  of  the  ions  is  low. 

Rather  than  introduce  new  parameters  we  shall  consider 
the  third  moments  to  be  0(e)  as  indicated  previously. 

We  now  take  higher  moments  of  the  ion  distribution 
function  and  express  these  in  Lagrangian  coordinates.   The 

computations  are  straightforward  and  will  be  omitted.   Upon 

2    2    2 
multiplication  of  (2.6)  by  u.^,  u-,^  \i    ,    \i    ,    u„,  u  ,  u-,u„j  and 
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u-,u^   we   obtain 
1   3 


l2 


r  v+^'    .    /^^+    N  T.   T.+ 


■%"] 


R"  +    (RPii)y    "    ^Y^22 


E   =    0(e) 


2^         2    ^    .^2     .    ^    ,^^+    N        .  1        ,^   5^2„^2, 


(r/)t-  $  \^    +  RY(RPri)T  + -2r-(V^  Pfi^Y 


(2.15)S 


T    ^T 

-    4   V%P+2 


^+ 


R\  Pii 


2RyR^PJ2-    2R^E 


0(e) 


_2.  /  o3.,+    ^+ 
R^ 


(V+2)t    +  I  RpV+^   +   RR^yP22    +   Ry^^P22)t   +  72^^^^""    Pr2^Y 

+   2R^RyP22   =    O(^) 
(¥+2)^   +    (RRyP^^)^   +   2(RW+P+2)y   =    O(^) 


2^+ 


(RV-^)^   +    (R^Pi2)y   =    O(^) 


¥+   +    (R   P+^)y    =    0(e) 


2.16)<^ 


(V+R^)^    +  ^  V+(R^)^    +   Ry(R   Pi2)t    +  -^-(R^Rt^Pi2) 


.+T 


,+ 


+T.       -r& 


+  (y-^R  p^^)y  +  v^Ry^ii 


EV 


+ 


R    R 
0(e) 


Y 


T 


(W%)^    +   Ry(RPi%)t    +  -2~--2(Rt  ^^i,)y 


r+T.+ 


+ 


r-r^p+^ 


13' 


+    (RW+PJ^)y   -    EW"^  =    0(e) 


The  right  hand  sides  of  (2.I5)  and  (2.16)  are  composed 
of  third  moments  and  terms  multiplied  by  e.   Notice  that 
the  first  equation  shows  that  E  is  indeed  of  order  one. 
The  four  equations  of  (2.l6)  show  that  if  V  ^  W  ,  P-j^2^ 
and  P^^  are  initially  zero  they  remain  0(e)  for  all  time, 
a  fact  which  we  already  knew  from  (2.1-^).   Discarding  these 
terms  from  the  left  hand  sides  of  (2.I5)  we  are  able  to 
eliminate  E  from  the  remaining  equations.   We  have  then 
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(RRy^pJi)^=  0(e)  ,     (r5RyPJ2)t  =    0{e)    , 
(2.17)  (RRyP^^)^^  =  0(e)  , 


'12  =  0(e)  ,         P+^  -  0(e)  ,    P+^  =  0(e)  . 


p. 


The  ion  pressure  terms  as  determined  by  (2.17) 
contain  three  arbitrary  functions  of  Y  which  must  be 
determined  from  initial  conditions.   Since  N  =  l/RRy 
we  can  write 

P+^  =  RVh^(Y)  ,   PJ2  =  -2  h2(Y)  ,   P35  -  N  h^(Y)  . 

¥e  shall  take  as  an  initial  condition  N  =  1  which  gives 
R  =  2Y.   Thus,  if  we  let  aQ(Y)  be  the  initial  value  of 
the  diagonal  terms  in  the  ion  pressure  tensor  we  find 
h-^(Y)  =  a"^/2Y,   h2(Y)  =  2a"^Y,  and  h  (y)  =  a"*"  leaving 

(2.18)      ^11  =  a"^/2YRRY  ,   PJ2  =  Sa'^Y/R^Ry  ,   P^^=  a'^ /RR^  . 

The  form  of  the  ion  pressure  tensor  is  unusual. 
P  R  P  R 

If  we  let  N*^=     rNdr/    r  dr,  that  is^  the  average 

number  density  of  the  Ions  between  the  center  and  R,  we 

can  express  (2.l8)  in  the  form  P^^  =  a''"N^/N*  ,  P22  =  a''"NN* , 

p";-,  =  a  N.   Near  the  center  N  will  not  differ  significantly 

from  Nj  assuming  the  fluctuations  In  density  are  not  too 

+       + 
violent.   Here  then,  both  P.-,  and  Ppp  are  nearly  proportional 

2  * 

to  N  .   Far  from  the  center,  however,  N  is  nearly  constant, 

depending  little  on  local  variations  in  density.   P, -,  now 

behaves  like  N-^  while  P„„  is  like  N.   R-^j  being  less 
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influenced  "by  the  cylindrical  geometry  is  everywhere 
proportional  to  N. 

c .   The  electron  pressure  tensor,  general  case 

To  determine  expressions  for  the  electron  pressure 
we  find  an  asymptotic  expansion  of  the  electron  distribu- 
tion function  and  from  this  determine  the  pressure  tensor. 
Our  method,  suggested  by  Gardner,  is  similar  to  that  used 
in  guiding  center  theory  but  differs  in  that  E,  is  0(l/£). 

Following  Morawetz  ( l4 )  we  write  the  electron  Vlasov 

equation  in  the  form 

g2  2 
fl   +  u^f;  +  ^  (— ""-2  -  BE^  -  su^B^  +  ^^562  )f-^ 

(2.19) 

1   ^^l^P  -      1 

2  3 

From  (2.19)  we  find  that  f~  will  be  constant  along  particle 
paths  given  by 

du-^  -j_   e  u^ 


dt 


=  U^    +  U^U^   =    —  {-J, BE^  -  BU^B.  +  BU^B^) 

Li  X  t. 


dUg  -^        BU-j^U^ 

dt"  =   ^2,  +  ^1^2   =  -  B  (--f-  +  ^2  -  ^1^5) 
t       r  ^ 


du^ 


-5-  -  u- 


dt    "3,     ^  ^ 

These  equations  suggest  looking  for  series  representations 
of  the  velocities  in  the  form 
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U  =  Uq  +  £U^  +  e^U^  +  . . . 

(2.21)  V"  =  I  V_^  +  Vq  +  s  V^  +  ... 
W"  =  —  V^_-|^  +  Wq  +  £  W-j^  +  .  .  . 

where,  motivated  by  (2.19)  we  define 

F  ^-1  "  ^-1^3  +  ^-A  -  E  =  0 

(2.22)  V_^   +  UqV_^   =  -  (-^--i  +  E^  -  UqB^) 

W_i   +  UqW_^   .  _  (e^  +  UqB^)  . 

We  define  new  velocity  variables  i,    r],    and  ^  by  |  =  u-,  -  Uq, 

T]  =  u„  -  —  V  .  ,  and  ^  =  ^^   -   -^  ^_1'   Clearly  ^  is  of  order 

one.   We  want  this  to  be  true  also  of  t]   and  C.   We  must 

therefore  restrict  the  values  of  u„  and  u^  for  which  f 

2      5 

is  non-zero  to  lie  almost  entirely  within  some  intervals 
of  V_-|  and  W  -,  .   These  intervals  may  be  large  but,  strictly 
speaking,  must  be  independent  of  e.   That  is,  the  electron 
temperature  must  be  uniformly  bounded. 

In  the  new  variables  (2.19)  takes  the  form 

b'f  df  ^^+  I  ^   l)^         1  Is 

(2.23)  -    B2(C+   ^  W_^)]f^    -    |[E2-(e+UQ)B^+  £   |(^+Uq)(ti+  |V_^)]^ 
-    I    [E^    -     (^+UQ)B2]f(;   =    0    . 

A  AAA/\ 

where  f  ( t,  r,  f„  n,  g  =  f  (t,r,u-^,U2,u  )  and  df/at=  f^+^^f^+^^f 
+  ^f^r   with  a  similar  formula  for   of/Or.   Using  (2.22) 
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this  can  be  written 

^x'^  >.'^        2V   ri  2 

(2.24)    +  [|  {V_-^   +  UqV-I   +  «Bj  -  I  ^)  +  a  (?+  Uo)]f^ 

L«  J. 

In  this  translated  velocity  space  we  define 
elliptical  coordinates  O^  9,    and  *  by 

•2  2         K.2  T^  ^  /o^fiN2  .  /^^    x2 


o2 


^  ^  +  -L^  +  ^5-  ,      tan   0   =  -J-  ,    tan   *  =  a_^^z. — ^--±i— 
p2        Q^        S^  ^^  (KiC)^ 


where  P,  Q  and  S  are  as  yet  undetermined  functions  of  r  and  t. 
O  is  a  measure  of  the  temperature  of  the  electrons.   These 
may  be  inverted  to  yield 

4=IOcos    0sin*j    ri^QPsin^sin*,    C   =   SOcos    4>. 

Letting  F{t  ,r,O,0  ,<i>)    =   t  {t  ,r ,  i,ri,  1^)    we  have  the  following 

identities '. 

r     -   cos    0   sin   tt>  sin   0      „  cos    0   cos    <i>   „ 

^i  ~  P  ^O   ~  QP  sin   CD  ^0   +  QP  ^D 

<^     _   sin   0   sin   4)   „  cos    0      „  sin   0   cos    <1)   „ 

^n    "  Q  >0   +  OQ    sin    0    ^0    "^  OQ  ^4) 

^  cos4)„  sin*^ 

^C   =  ~^        ^O    "   T)S        ^* 

Sf  tP?  '^t??  tP 

^-  =   F^    -   F^[0  —  cos    0    sin    4)    +  O  ^  sin   0    sin    *   +  O  -^  cos^* 

^t  "^t  ^t 

-   ^5-  cos    0   sin   4)   -  -p--  sin   0   sin   4)    -  -g-  cos    *  ] 

P        Q,  ?,sin   0        T]   cos    0 

+    F0t(p--Q-)^i^    0    COS     0     -    ^-^-^   +   ^\-^^,]  + 
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+  F^[f^  Ucos  e    (SQ^)^  +  sin  6    (SPri)^)Gos  (f 

-  (PQC)^  sin  <t>}  ]  .  , 

A  similar  formula  holds  for  of /or  . 

When  these  differential  identities  are  substituted 
into  (2.24)  we  find  that  by  setting 


2  rB^  -  2V  ,  p2        B. 


the  coefficient  of  —  F^  vanishes.   Equation  (2.24)  becomes 

F^  +  (P^  cos  9  sin  *  +  Uq)F^  +  ^  F^  +    ij   €  '    +jC)')^e 
(2.26) 

where 

2   Q 
^  =  -  Q^P^cos^e  sin^O  +  O^P(p  ^  -  j-^  -  ^)sin^0  cos  9  sin^* 

-  0(/-  +  Uq  /  +  UQ^)cos^e  sin^*  -  0(q1  +  Uq  ^  +  ^) 

P  p  ^-j-  r2  2r  2 

•sin   0   sin   <t>   -  n(^-  +  U^  ^— )cos    *    -  O     — <^  cos    9   sin   4>  cos  <1> 

-  p(Uq     +  UqUq    )    cos    0   sin   * 

t  r 

rj  ,         F    /    -n           I    ^T      \    \        P    /-D      ,    Tf        \    sin   0   cos    * 
e'    =   Q    (rB^    -    (rV_^)^)    -    s     (Bp    +  W_^^) ^^-^-^ 

c^'    =  —    (Un   +  UnU^    )    -!"    g   +    (-^  -   ^  +  U^(-^  -   ^t) 
OP        ^t         0   Op      sm    4)  P  Q  O'P  Q    ' 

^  +  U^    )    sin   0   cos    0   +  OP(--i-  -  ^r -  -  -)sin  0  cos''0   sin   <t> 

r  0    '  ^P  Q,  r' 

2 

-  -^-  sin''©    sin   4) 

rP 

£-   =      I    (B^    +  W_,    ) 

r 
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S,    P,       S     P 
^"  =  (g^  -  /  +  Uq(3^  -  /))cos  *  sin  * 

+  OP(7=^  +  p— •)  cos-^Q  sln-^*  cos  * 

-  QS(-p-  +  Q— )  cos  6   sin  <t)  cos  * 

QPp   QSp    0^    1      2  2 

+  OP( — o  +  CQ—  +  ~^~  ~   ^)  sin  9  cos  0  sin  4>  cos  * 
P^    ^^     P^r    ^ 

1  2 

-  ^=-^   {IJ        +  ^JJq    )    COS    Q   cos    *    -   Upj  cos    0   cos    <t>   sin   <t> 

Uq  2 

-  —  sin    0    cos    0    sin    4> . 

r 

Equation  (2.26)  appears  to  have  several  singularities 
occurring  at  r  =  0.   However^  if  we  take  the  natural  condi- 
tions that  V~  and  U  vanish  at  the  origin  we  find  that  the 

1  0^ 
term  —  (-^^  -  l)  which  appears  in  ^   and  /S)'    can  be  written  as 

r  (p?  -  1)  -  r  (rBj-TTV:-r;  -  ^'  -  Bj"+  t(rv:j;   ' 

This  will  be  regular  if  P/Q,  is.   We  shall  later  solve 
for  P,  Q  and  S  explicitly  and  verify  that  this  quotient 
is  well  behaved.   Since  U^^/r  is  regular  the  only  true 
singularity  appears  in  ^'  giving  the  condition  that 

Fg(t, 0,0,0,*)  =  0. 

We  seek  now  a  series  representation  of  F  in  the  form 
F  =  Fq  +  e  F-^  +  e^F^  +  .  .  .  . 

This  is  substituted  into  (2.26)  and  the  coefficient  of 
each  power  of  e  set  equal  to  zero.  The  first  equation 
yielded  is 
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(2.27)      i'sir.   5  ^^  ^  r"  [ZZB    i  sin  *  i^Q  -   sin  a  cos    #■  F      }   =    D 
nZae   ejLxii^i^    siluTizr.   of    ^2.27)   is   ?-   =    7  z.r.n,^')    li'neyre 


-f  1  -      ■         -  1  --     \ 


distribution  i^anatlcn  Is   IndeDeniienT    cr    Z  =;    -nat 


-,n)    . 


^_      ^   (?n  COS   ?   sin   ff   -5-^0)^     ^   ^^=^_    ~   ~'^- 
(2.25) 

-=-  r"[sos   5  F,       -    ':~    \    £ln   r   J-     '    =    C 
9  - 

iSils   1=   Inxegre'ei   »ltn  respect  to  5  frcci  r  =   0. 

t  T 


-  ^n   )    t^  -  ^--—      ~^~-^^ 


0^ 


(2.29) 


13C 


S^  s 


-F  r.~    — .^^    sin    r    sin    ?     2DS    f    ^  '^{-^  4-    Up,   ^s-^)5    COS"*    4 


_c^_ 


+  i   (Uq      +  UqUq    )    sin   9   sin   <t>]^Q 
t  r      g  O 


=    _  r-F,   -  r"  (4-  r  cos  e  F.  ae  +  ^^^  T  cos  ©f, 

1  ^o<i>  J  1  sxn*^  1 


de 


0 


cos    4) 
sin 


I  sin   0F^)    +    ^-^{t,r,n,'i>) 


We  make   two   definitions  I 

27r  2ir 

L(t,r,0,*)    =   J  cos    0   F-j_   dQ,  L'  (t,  r,0,  *)    =  J  cos    9    sin   0F^d0, 

0  0 

Using  the   fact    that   F-,    must  be  periodic   in  9  we   evaluate 

(2.29)  at    0  =   2nTr^    where   n   is   an   integer  to   obtain 

V  ^  Pf        ^t  ^r        "^r 

2nTr(5^     +  Uq  /q    )    -  OnTr[(/  +  _^-)    +  U^l/-  +  ^^-) 

t  r  ^ 

U  S  S 

(2.30)  +  Uq   +  -^)sin^4)   +   2(g-^  +  U     -^)cos^(t)]S^Q 

r  O 

=      -    r.F^(t,r,0,2nTr,*)    _  nP' (l^+  ff-^-^  L)+ i?-^  (t,  r,  *) . 

Setting  n   =    0  gives    5^-,  ( t,  r^O,  *)    =   F-,  (t ^  r,Q,  0^  *) .      Then 
for  n  nonzero  we   obtain   from  equation    (2.30) 

2(5^0/  Uq  5^0    )    -  0[(p^  +  ^  +  Uo(/  +  Q^)    +  Uq   +  -^sm^c 

0  c 

+  2(3^  +  Uq  /-)    cos2,]5^o^   .    .   ^L._  |_   (L   Sin   *). 


2.31) 


Multiplying  this  by  sin  <t>  and  integrating  with  respect  to 
from  <t)  =  0  gives 
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sin  *L  =  p-  cos  4>  ( SC  +  Uq  5^  ) 

t       r 


+  pr  [(p-  +  ^  +  Uq(^  +  ^)  +  Uq  +  — )(-^ cos  *) 

(2.32)   -  2(k-  +  Uq  S^)  --|-]>c^-  g(t.r,0)   where 


^r   ^r   ^r         ^0  c^ 

Evaluating  (2.52)  at  O  =  tt  gives  simply  2g  =  0,  so  that^ 
since  g  is  independent  of  <t>,  it  is  identically  zero.   We 
still  have  one  relationship  in  P^  Q  and  S  at  our  disposal. 
By  taking 

(2.33)   PQS  =  A(Y)N,  with  A  an  undetermined  function  of  Y;, 

the  second  term  in  g  vanishes  leaving 

(2.3^)  5^0  +^0^0  =  °  • 

t        r 

Equation  (2.32)  with  g  set  equal  to  zero  gives  an 

explicit  expression  for  L.   This  will  be  needed  to  determine 

moments  of  the  distribution  function.   We  shall  also  need  an 

expression  for  L' .   For  this  we  multiply  (2.28)  by  sin  9 

and  integrate  with  respect  to  0   from  0  to  27r^  taking  (2.32) 

into  account.   After  integrating  by  parts  we  obtain 

(2.35)  r"(sin^<fL' )^  =      T'    sin   <t>L    . 

Upon   solving  for  L  and  L'    we   obtain 
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(2.36) 


From  equations  (2.25)  and  (2.33)  we  find 
g       ^  rB-,  -  2V  ^ 


^'  -  (^^)'  FB-|^rrx-iT 


•l^r   ^2+^.1 


.  rB  -(rV  .)   -2 

(2.57)  Q^  =  (AN) 2  (^^^TW— -) 


B. 


rB^  -  2V_^ 


B^  +W_-j^ 


B. 


^  rB   -(rV   ) 

P  r 


•2 


We  are  now  able  to  evaluate^  to  lowest  order  in  s, 
the  moments  of  the  electron  distribution  function.   For 


the  mean  radial  velocity  we  find 

00  27r  TT 


U-,  f      I  du  j 


J  J  J^"^^  ^^^'^  ^^^  -tJ+U  JO^PQS    sin^OF  d*   dO  d( 


Jf-  Ida 


OOP 
CX)  27r  TT 

in 

000 


O^PQS    sin^tDFd*   d0   dO 


00     2Tr   TT 

'0    '    N  J       J  ^ 

0      0      0 


=   U^   +  ^    I       I       I   O^   COS    0   sin^*    (5^^p)+   eF.  )    d*   d0   dO  +  O(e^) 


.3. 


^0 

=  Uq  +  -w- 


CO      21^      IT 


J     J     O^   cos    0   sin^4)   d*   d0   dO 
0      0      0 


00      TT 


e 
+  N     .. 


C?   sin^ct)   L  d*   dO  +  0(£^) 


0      0 
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The   first  integral  can  be  evaluated  directly  and  is  seen 
to  vanish.   Using  the  expression  for  L  in  (2.36)  the  second 
integral  is  seen  also  to  vanish.   The  remaining  velocities 
are  found  in  a  similar  manner  yielding 


1 


1 


(2.38)   U  =.  Uq  +  O(e^),   V"  =  i  V_^  +  0(e)  ,   W"  =  ^  W_^+  0(e) 

Using    (2.38)    the   electron  pressure   terms   needed   in 
the   equations    of  motion   are   evaluated! 
2 


Pll  =  -^  J  (u,-u)2f-  |da| 


^1" 


00     27r     TT 


l  +  £ 


-2 


[P%  J  J 


O^cos^   e   sin   4>F  d*   dS   dQ  +   0(e^)] 


0      0      0 


Pop  ^  ---p    ^   (up-V 


du 


12 


l+£ 

l  +  £ 
2 


00      2F      TT 

2    [PQ^S  J     J     J  O^^   sin^e   sin^(t)Fd*   d9   dO   +  0(e)] 
0      0      0 


2 


lu 


P    [P  Q  o 

1  +  e^ 


,2 


(u^-U)(u2-V    )f      Idi 

00  27r  TT 

'  J  J  ^  "^^^    ^   ^^^   ^   sin  %  F  d*  d9  dO  +  0(  e^)  ] 
0   0   0 


— ^    r   (u. -U)(u^-W    )f      |du| 
P  00  2Tr  n- 


=      -~o    [P  QS  O      COS    0   COS    <J>    sin-^(t)F  d*   d9   dO+0(£    )] 

Ire  J   J   J 

"^^  0  0  0 

When  the  lowest  order  expression  for  F  is  substituted  into 
these  expressions  we  can  integrate  to  obtain 
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CO 


11 


22 


2  ll_  p^Qs  r  a^ 5^   dO  +  o(£^) 


"5" 


,2  57r 


0 


0 
oo 


"8" 


PQ-'S 


OD  TT 


0 


^22 
P;L2  =  £  P  Q  s 


O^l'  sin  *  d*  dO  +  0(e) 


(2.59) 


0  0 


00 


5  1^4'  (s.+usJpV  ro^5^,  dO+0(s^ 


9r' 


,2   ^^t 

00  TT 

.4 


0 


P"   =  E^P^QS^  J  J^  L  cos  4)  sin^*  d*  dO  +  O(e^) 


13 


0  0 
.2 


CO 


=  -  e 


3  _5Ti:  (s  +uS^)p2QS  f  O^  5>.  dO  +  0(£^ 
l6r"   t     r      0 


,2 


where  we  have  used  (2.36)  a^^d  (2.38). 

It  is  necessary  here  to  discuss  the  relative 
magnitudes  of  the  ion  and  electron  pressure  terms.   The  e' 
multiplying  the  electron  terms  would  seem  to  indicate  that 
these  terms  are  small  compared  to  the  ion  terms  and  that 
since  we  have  already  neglected  terms  of  order  e  we  should 
neglect  these  as  well.   This  is^  however,  not  realistic. 
The  average  kinetic  energy  per  particle  must  be  the  same 
for  the  electrons  and  ions.   In  a  collisional  theory  this 
would  be  true  and  although  our  model  is  collisionless  it 
can  be  thought  of  as  arising  by  allowing  a  collision 
parameter  to  go  to  zero.   To  be  consistent  with  our  asymptotic 
development  and  yet  retain  both  pressures  in  the  equations  it 


,00 


is  necessary  to  assume  that  the  quantity  £     O  J- ^   dO  is, 

2        ^0     ^ 
although  strictly  speaking  stil  0(e  ),  sufficiently  large 
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to  allow  the  electron  pressure  terms  to  be  comparable  with 

p^  2  ^ 
the  ion  terms.   Since    Q  :f^   dD  is  the  number  density^ 

and  thus  approximately  one,   we  must  assume  that  the 

thermal  speed  of  the  electrons  is  high.   We  define  the 

quantity  a"  to  be 

(2.40)         a"(Y)  =  e^  ^-  J  o'^  3^^  dO  , 

and  assume  that  a"  is  sufficiently  large  that  it  cannot  be 
neglected  compared  to  one.   With  this  assumption  the 
electron  pressure  terms  are  of  the  same  magnitude  as  the 
ion  terms.   Equation  (2.34)  shows  that  a"  is  indeed 
independent  of  T.   We  then  have 

^11  =  ""P"^^  S  +  O(e^) 

P22  =  ct"P  Q^S  +  O(e^) 

(2.41) 

P"   =  -  £  ^—%   (S,  +  US^)  a-pV  +  O(e^) 
12       27Trr""^    ^      ^ 

P"   =  -  e  -^,7  (S,  +  U  S  )  a'P^Q  S  +  0(£'^)  . 

d.   The  electron  pressure  tensor,  theta  and  pure  z-pinches 

To  study  the  theta  pinch  in  which  'Q   ,    ^   ,    W  ,    and  W 

are  all  identically  zero  we  observe  that  since  there  is  now 

no  force  in  the  z  direction,  S  is  identically  one.   By 

setting  S  =  1  in  (2.57)  we  see  that  in  the  limit  where 

Bp  and  W~  tend  to  zero  we  have 

Bp      2  rB   -  2  V  . 

(b--4:~)  —  (^^^  rB-!---c7v:;T;  • 

r 
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Using  this  we  obtain  from  (2.37) 

(2.42) 

u  „   rB^  -  2  V  .   -1 

5     -1 '  r 


P~   and  P~   are  now  easily  found  from  (2.42)  and 
(2.4l)  with  S  set  equal  to  one.   To  find  P~   we  must  refer 
back  to  (2.28)  which  now  becomes 

(2.43)  5C  +  (PO  cos  0  sin  *  +  u^)  ^j'^   +  -^^5^^  +  ^'F.   =  0  . 

^t  r       Q      0 

Using  the  fact  that  (B^/i\l)r|^  =  0  we  integrate  (2.4^)  with 
respect  to  0  and  find  F-,  .   From  this  the  pressure  terms 
appropriate  for  the  theta  pinch  are  found  to  be 

P~^  =  a"P^Q  +  0(6^) 

(2.44)  P22  =  a"P  Q^+  0(6^) 

P12  -  ^  I5-  -^~-   «"  (^  +  U  ^  +  U^)  +  0(e^  . 

P~   is  no  longer  needed  in  the  equations  of  motion. 

When  B^j  E„,  V"*",  and  V~  are  all  identically  zero 

j?    2 

we  have  a  pure  z  pinch.   By  setting  Q  equal  to  one  we  find 

that  in  the  limit  when  B^  goes  to  zero 

rB^  -  2  V  T    2  ^     B„ 

^B^-  (^V_3_)^  B2+W_3_^ 

For  this  case  we  obtain  from  (2.37) 
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P^  =  (AN)^      ^ 


(2.45)  '  "" 

r 
As  before,  we  find  P"   and  P~   from  (2.45)  and  (2.4l) 
with  now  Q  set  equal  to  one.   P-,  ^  requires  no  special  care 
in  this  limit.   The  electron  pressure  terms  for  the  pure 
z  pinch  are  given  by 


P"   =  a  P^S   +  O(e^) 


11 


(2.46)  P^^  =  a  P  S   +  O(e^) 


22 
Pl3  -  -  £  I   -^^  a"  (S^  +  U  3^)  +  0(£^) 

For  the  pure  z  pinch  P-,  p  is  not  needed. 

e .   Th e  s y s t ems to  be  solved 

Having  found  expressions  for  the  ion  and  electron 
pressure  tensors  the  equations  of  motion  (2.15)  and  (2.1^^) 
are  now  closed.   To  incorporate  the  electron  pressure  terms 
into  the  equations  of  motion  we  must  first  express  them  in 
Lagrangian  coordinates.   We  then  find 


H 


P  ^  -^^2  ~  (^^2^YY 


H^  -=    m^B^    -    (R^B^Y^Y 

,    ,  2 

P^  .  (AN)^  (1  +  — ^")  (1  +  — /--)  +  0(e) 

"2  "3 

(2.^17) 


2 
S^  =  (AN)2  (1  +  --4^)'  (1  +  ~r^^-)    +   0(e) 

o 
Q^  =  (AN)2  (1  + 1-^^)  (1  +  --TT^-)'   +  0(e) 

T'  =  P^S  E-,    -  0(e) 

A^  =  0 

When  these  terais  and  the  expressions  for  the  ion  pressure 
are  substituted  into  (2.13)  and  (2.l4)  we  obtain,  after 
eliminating  Q,  the  following  closed  system  of  equations. 


I^,^    -    R(B^y)^    +  B^iBB^)^   +  R   B^B^Y 


_^    /    a+       X         2a  T    ,     ,a  AP^s  a  AQ"  ,    . 


(2.48)  2YRy-Y        j,3  Ry      Y  R 

[T       a   ASm 
\    -   T2    (-H--)yY=    °^^) 

For    the    theta    pinch   in   which   only   the   B^   field    is 

present   we    obtain 

2 

p^    =      {mf  (1  +  -4"-)  +  o(^-) 

2 
Q^    =  (AN)2    (1    +  — H^-)         +   0(e) 

T'    =  P^H      +   0(e) 


(2.49) 


A^  -  0 
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The  equations  of  motion  become 

-   2       -   2     "^ 

a  AP      a  AQ 

(2.50)  +  ^-R7-)y  -  --R-  =   °(^) 

_  2 

And  finally,  the  pure  z  pinch  expressions  become 

P^   =   (AN)^  (1  +  --^^^-)+  0(£) 

"2 

S^*   =   (AN)^  (1  +i^^Y-)'^  +  0(e) 

(2.51)  pH  ^ 

A^   -   0   , 
and  the  corresponding  equations  of  motion  are 

(2.52)  +  (24^-)y  -'^     =   0(=) 

2 

2^    12     ^2   YY 
Systems  {2j\Q),    (2.5O)  and  (2.52)  are  solved  with 
the  following  initial  and  boundary  conditions.   We  assume 
that  initially  there  is  an  infinite  cylinder  of  plasma  at 
equilibrium  in  which  the  density  is  everywhere  equal  to  one 
and  there  are  no  currents  in  the  plasma   Thus,  B^  is 
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identically  zero  while  B^  is  constant.   To  initiate  the 

motion  an  increasing  magnetic  field  is  applied  at  the 

outer  boundary  of  the  plasma  causing  the  radial  contraction. 

When  pressure  is  added  to  the  problem^  however^  we 
find  that  we  have  more  free  variables  than  vje  need,  for 
we  now  have  a(Y),  and  A(Y)  which  are  as  yet  undetermined. 
If  we  assume  that  the  initial  electron  distribution  function 
is  symmetric  about  the  mean  velocity  we  have  by  definition 
that  P,  Q  and  S  are  constant.   From  (2.33)  we  then  find 
that  A(y)  is  also  constant.   Without  loss  of  generality 
we  take  these  constants  to  be  one.   When  the  desired  initial 
values  of  R,  Bp  and  B^  are  substituted  into  the  equilibrium 
equation  obtained  from  setting  all  T  derivatives  equal  to 
zero  in  (2,48)  we  find  that 

2Y(a"^  +  a~)^   -  0  . 

Since  these  terms  are  positive,  they  must  be  constant. 

The  initial  and  boundary  conditions  for  the  z  and 
theta  pinches  associated  with  systems  (2.48),  (2.5O),  and 
(2.52)  are  discussed  in  Appendix  5.   They  are 

N  =  ^  =  1  =>  R(0,Y)  =  v^Y 

R^(-0,Y)  =  B2^(-0,Y)  =  B^^(-0,Y)  =  0 

65(0, Y)  =  0 

fl   for  (2.48)  and  (2.50) 
B.(0,Y)  = 


^  0   for  (2.52) 

(2.55) 

R(T,0)   =   0 
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B    (TjO)      regular 

B2(T,Y^)    prescribed      ;  (RPll)v    "   ^Y^22 


B^(T,Y^)   or      B^    (T,Y^)  prescribed. 


=    0 
Y=Y 


Y   is  the  value  of  Y  at  the  outer  boundary  of  the  plasma. 
m 

We  shall  consider  also  some  pressure-free  flows  in 

which  both  fields  are  initially  present  in  a  force-free 

1   2 
configuration.   Let  J^(0,Y)  be  a  constant  denoted  by  ^^  ^^  • 

With  this  definition  B2(0,Y)  -  mY  and  (2.53)  must  be  altered 
by  letting 


B^lO.Y)  =  mY   and 


(2.53a)  B^(0,Y)  satisfy  (B^  ^)      -   B^B^   -  m  -  0 
-^ith  B^(0,Y^)  +  B^(0,Y^)  =  1  . 

The  quantities  Hp  and  H^  defined  in  (2.47)  are  seen 

to  arise  in  two  different  contexts^  indicating  that  they 

have  significance  in  their  own  right.   This  notion  is 

enforced  by  the  fact  that  when  the  electron  pressure  is 

absent  they  are  both  invariant  in  time.   Their  relationship 

to  tlie  more  familiar  one-fluid  fields  can  be  seen  by  writing 

them  as 

B  B 

^2  =  W  -  ^^  ^2)yY  ^""^      ^3  ^  N     ^^  V^Y  • 

Thus,  when  the  differentiated  terms  on  the  right  are  small, 
or  when  they  vanish  as  they  do  initially,  the  near  invariance 
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of  these  terms  says  that  Bp  <X.  RN  and  B^cCN.   These  are  the 
expressions  for  the  fields  in  the  one-fluid  model.   From 
(2.48)  we  see  also  that  the  two  field  equations  are  coupled 
by  these  quantities. 

f .   Characteristics 

It  is  clear  that  the  inclusion  of  electron  pressure 
greatly  complicates  the  equations  of  motion.   It  also  raises 
some  new  questions.   For  instance,  since  the  highest  deriva- 
tives now  appear  with  a  factor  of  a  do  the  equations  become 
singular  as  this  parameter  goes   to  zero?  Furthermore,  from 
the  definitions  of  P  given  in  (2.4/),  (2.49),  and  (2.51) 
for  the  different  variations  of  the  problem  we  see  that 
while  in  (2.48)  and  (2.52)  the  first  equation  contains  R^^ 
and  Kryy   as  the  highest  derivatives,  the  corresponding  equa- 
tion in  (2.50)  contains  only  Rrnm  and  Ryy-   Does  this  mean 
that  the  two  systems  are  of  different  type?  Finally,  are 
the  three  systems  with  initial  and  boundary  conditions  as 
given  in  (2.53)  and  (2.55a)  well  posed?  To  answer  these 
questions  we  must  examine  the  characteristics  of  the  three 
systems . 

The  characteristics  of  the  system  (2.48)  and  the 
special  cases  derived  from  it  are  determined  in  Appendix  1. 
The  results  are  given  here.   Let 

C+  =  ^-^  and    C-   =  ^^-^   . 
^     2YRy  ^       Ry 


■59- 


These  are  the  sound  speeds  corresponding  to  the  two 
pressures.   They  are  given  in  the  Lagrangian  representation, 
dY/dT.   To  obtain  the  more  familiar  Eulerian  form  it  is 
necessary  to  multiply  each  by  R^.   The  characteristics 
for  the  three  systems  with  and  without  pressure  are  given  by: 

Pure  z  pinch -system  (2.52) 


ct"^  -  0 
a"  =  0 


2    ^ 

dT   dY^  =  0 


a^  >  0 


=  0 


dT^  dY   [dY^  -    C'^     dT^  ]  =  0 


a"*"  >  Ol 


a 


>oI 


a 


,+ 


0 


a  =  0 


4        2 
dT  dY   [dY^^ 


.+ 
^L 


c;   +  C^  )  dT  ]  =  0 


Theta  pinch  -  system  (2.50) 


2    ^ 
dT   dY^ 


0 


a+  >  0 


a 


0 


dT^  dY   [dY^ 


cj  dT^]  =  0 


a"^  >  0 


a"  >  0 


4        2 
dT   dY   [dY 


n  + 


.2 


(Cj   +  G^  )  dT^]  =  0 


Both  fields  -   system  (2.48) 


a"*"  =  0 


a"  =  0 

+ 


>  0 


=  0 


4   4 
dT   dY 


0 


^1      P       P       4-       ? 

dT'  dY   [dY^^  -  Cj  dT""]  =  0 


a"^  >  0 


a   >  0 


62     p      <^  _^     p 

dT   dY   [dY^^  -  ('^7   +  C.  )  dT^] 


L 


L 


0 
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By  the  nature  of  the  characteristics  we  see  that 
the  prolDlems  are  essentially  fluid-dynamical  in  nature, 
the  only  non-trivial  characteristics  being  the  sound  speeds. 
The  electron  sound  speed,  like  the  electron  pressure, 
depends  on  the  mag.netic  fields  as  well  as  on  the  density. 
Notice  that  the  Alfven  speed  does  not  occur  as  a  character- 
istic as  it  does  in  the  one-fluid  magnetohydrodynamic  model. 
¥e  see  also  that  the  characteristic  equations  are  -well 
behaved  when  the  electron  pressure  goes  to  zero.   The 
vanishing  of  a"  causes  only  the  disappearance  of  the  electron 
sound  speed  from  the  characteristic  equations.   Also,  as  was 
to  be  expected,  there  is  no  basic  difference  in  the  character- 
istics between  the  theta  pinch  and  the  other  two  problems. 
The  third  derivatives  in  R  which  appear  in  the  momentum 
equation  of  these  systems  is  eliminated  by  virtue  of  the 
fact  that  the  highest  order  parts  of  the  field  equations 
are  not  independent  of  these  terms.   Finally,  since  the 
systems  have  only  these  fluid-dynamical  characteristics, 
the  systems  (2.48),  (2.5O)  and  (2.5?)  are  well  posed  with 
the  initial  and  boundary  conditions  given  in  (2.53)  and 
(2.55a). 

g ,   An  energy  integral 

Wlien  the  electron  pressure  is  absent  and  a   is  taken 
to  be  constant  an  energy  integral  is  associated  with  the 
system  (2.48).   By  combining  the  equations  in  the  following 
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manner 


2R^  [R^-^    -    R{B^^f   +   RB^B^Y  +  B2(RB2)y 

+  a+(— -^)^  -  ^^^}  -  2B   {(rS   )   -  RR^B  } 
2YR^        R^       -^       -^         i  ;?  ^ 


-  2R32  {(RB2)yy  -  Ry^2^T  "  ° 


we  arrive  at 

[r^  +  rRy(b^  +  b|)  +  rS^y^+  (RB2)y 

+      + 
(2.5/1)    +  _«_^  +  -^}   +  {RR^(B^  +  b|)  -  2B^(R^B^y)t 
2YRy    R   T 

"P 

-  2RB2(RB2)yt  +  a"^  ~^^Y  "  °  ' 

yRy 

This  integral  provides  a  useful  check  on  the  accuracy  of 

the  numerical  calculations  and  serves  as  a  guide  in  the  proof 

of  the  stability  of  the  finite  difference  scheme. 


h .   The  one -fluid  model 

In  this  section  we  compare  our  two-fluid  model  of 
the  pinch  with  the  simpler  one-fluid  model.  To  do  this 
we  write  the  Lundquist  equations  in  axially  symmetric 
cylindrical  Lagrangian  coordinates.   Let  N  he    the  total 
fluid  number  density,  V  and  W  the  circumferential  and  axial 
velocities.   The  pressure  is  taken  to  be  proportional  to 
the  density  squared.   Keeping  the  same  scaling  as  in  (2.4) 
the  equations  of  motion  are 
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(Ry-^2^T 


0 


(RRyB  )   =  0 

(2.55)  ^ 

(RV)^  =  0 
W^  =  0 
(p(rRy)^)^  =  0  . 

In  comparing  the  two  models  we  shall  take  the 
pressure  to  be  zero  in  both.   In  the  two-fluid  model  this 
requires  that  the  net  momentum  in  the  9   and  z  directions 
remains  identically  zero^  the  value  it  had  initially.   To 
correspond  to  this  we  take  V  and  W  to  be  zero.   To  facili- 
tate comparison  we  rescale  the  independeiit  variables  in  both 
models  by  defining 

X  =  5^Y  ,  T  =  5  T 

where  5  is  a  positive  number.   The  two  fluid  equations 
of  motion  as  given  in  (2.48)  become 

(2.56)  (RR^B^),  =  s'^rSx^xt 

(R^B^)^   =   5^(RB2)XXT 
while  the  one-fluid  equations  (2.55)  are 
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(2.57)  (RRxB^^T  =  ° 

Upon  examining  these  two  systems  it  is  clear  that  as  5 
approaches  zero  the  two  fluid  system  (2.56)  goes  over 
into  (2.57)^  the  one-fluid  system.   The  limiting  process 
iSj  however,  singular  as  6  multiplies  the  highest  deriva- 
tives in  the  two-fluid  field  equations.   Thus,  as  Morton 
(15)  points  out,  the  two-fluid  model  does  not  allow 
discontinuities  in  the  mag.netic  fields  while  the  one-fluid 
does.   Both  allow  jumps  in  the  fluid  variables  when 
pressure  is  included.   As  6  goes  to  zero  any  oscillations 
present  in  the  two-fluid  model  become  compressed  both  in 
space  and  time  so  the  limit  may  be  though  of  as  a  space- 
time  average  of  the  oscillation. 

Continuing  the  correspondence  between  the  two  models 
we  again  take  the  initial  condition  to  be  one  of  equilibrium. 
In  this  case,  though,  we  can  find  simple  formulas  for  the 
initial  fields.   Thus 

(2.58)  R^(Y,0)  =  2Y  ,  Ti^{YO)    =   m^Y  ,   B^(YO)  =  b-  2m2Y  , 
where  m  and  b  are  constants.   From  (2.55)  and  (2.58)  we  have 


R   _  1   m     ani     R   ^  /b-2m2Y 


_i|2|- 


From  these  we  can  find  one  second  order  equation  for  r: 

(2.59)   RtT  +  I  «  t^f^-|]Y  +  ^   (/)y  +  AR  [— ^-plv  =  0 

The  characteristics  of  system  (2.55)  are  given  by 

d^^  [{^f   -  ^  (B^  +  b|  +  2P)]  =  0  . 

The  last  term  in  this  is  the  sum  of  the  squares  of  the 
Alfven  velocity  and  the  sound  speed  in  Lagrangian 
representation. 
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CHAPTER  III.   NUMERICAL  METHODS 

The  solutions  of  the  systems  (2.48)  and  (2.5O)  are 
studied  numerically.   Since  these  systems  are  of  mixed 
hyperbolic-parabolic  type  we  expect  the  stability 
criterion  of  the  difference  scheme  to  depend  not  only  on 
the  sound  speed  but  on  other  parameters  as  well.   The 
complexity  of  the  systems  makes  it  doubtful  that  such  a 
criterion  can  be  found.   It  is  possible^  however,  to  find 
a  stability  criterion  for  the  simpler  system  in  which  the 
electron  pressure  terms  are  absent.   The  existence  of  the 
energy  integral  allows  us  to  apply  the  so  called  energy 
method  to  this  case.   Having  done  this  we  then  solve  the 
full  system  and  compare  the  results  with  those  known  to 
be  stable.   If  the  addition  of  a  small  amount  of  electron 
pressure  does  not  cause  seriouo  disturbances  in  the  solution 
we  can  assume  that  stability  has  not  been  adversely  affected 
by  the  addition. 

a .   The  solution  without  electron  pressure 

To  prove  the  stability  of  the  finite  difference 
equations  used  to  solve  (2.48)  with  a   equal  zero  and  a 
constant  we  use  the  energy  method.   This  entails  finding 
a  discrete  analog  of  the  energy  integral  (2.54).   By 
determining  when  the  growth  of  this  quantity  does  not 
exceed  that  allowed  by  the  accuracy  of  the  scheme  we  arrive 
at  a  criterion  for  stability.   To  obtain  this  ^liscrete 
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energy  integral  it  is  necessary  to  make  some  preliminary 
changes  in  the  equations  of  motion  so  that  we  are  actually 
solving  a  slightly  altered,  though  equivalent,  system.   We 
first  scale  the  independent  variable  Y  by  a  factor  A.   This 
allows  us  to  eadily  change  the  initial  radius  of  the  plasma 
column.   (2.48)  is  then  written  in  the  form 

A  R 


(5.1) 


^^%yVt  =  I  [(R^)y^3]t 

R(RB2)YyT  +  I  R"  ^^^)y^2  "  "I  (R^)yT^2  "  I  ^^^)y^2T  =  ° 


R. 


i^v  ^  'O      (f^vZ)rn  +   (RmZ), 


The  third  equation  results  from  multiplying  the  second  member 
of  (2.48)  by  R  and  using  some  elementary  differential 
identities.   As  is  clear  by  comparing  (3.1)  with  (2.48),  Z 
equals  1/(YRy  ).  We  represent  Z  by  the  last  equation  rather 
than  by  its  more  logical  definition  because  it  will  be  our 
aim  to  express  the  change  in  the  energy  integral  in  terns  of 
the  T  increment  and  the  explicit  appearance  of  Y  makes  this 
im.possible.   This  will  be  discussed  in  detail  later. 

To  approximate  (3.1)  with  difference  equations  the 
(T,Y)  plane  is  divided  into  a  rectangular  lattice  with  sides 
parallel  to  the  axes  and  dividing  the  T  and  Y  axes  into 
intervals  of  length  k  and  h  respectively.   The  axes  form 
part  of  the  lattice.   To  minimize  confusion  with  subscripts 
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let  B  represent  B,  and  M  represent  Bp.   Define  the  differ- 
ence variables  r'^,  m'?  and  Z^  as  R(nkjjh)j  M(nk,jh),  and 
Z(nk,jh)  while  B^!  will  be  B(nk,  ( j-l/2)h) .   We  will  also  have 

J 

to  define  some  difference  operators! 

A^  -  a"}-^               ^     A^+^  -  A^       ^     A^+^  -A^^-l 
5-An  _  __J_ 3 5+.^  _  _3_ 3_  gO^n   _j_ j 

A  A^  -      j^^      ,    A  Aj  -      ^       ,    A  A^.  -      j^ 

oA.  =  ooA.  ,    AA^=AAA.  . 
J        J        J         J 

The  numerical  solution  of  (5.1)  proceeds  as  follows. 
From  (2.53)  and  (2.53a)  we  obtain 


J  -  -J  -  ^-3       '         3    -      3    -        ^-^  '         J  -  -J  -  jh(AQRl)5 

2 
(3.2)   B^  =  B^   where   (A+B^)^  -  A"'"B^   -  Am^  =  0 

2 

R^  =  M^  =  Z2  =  0  ,    bJ   =  g^(nk)  ,    M^   =  g^  (nk),A°R^=0, 


for  1  Jl  j  ^  J.   J  is  the  number  of  mesh  points  along  the  Y 
axis  and  is  taken  to  be  400  in  our  problem.   Assume  now  that 
we  have  R^  M^  B  and  Z  for  all  j  at  n-1  and  n.   (3.2)  gives 
this  for  n  =  2.   r"^    is  found  from  the  first  equation  of 

J 

(3.1)    discretized   in   the   form 

R^+1    =    2R^    -    R^-1    +  k2{    1      R^(aX)^    -    A  R?(A+Bf  +   ^^^^^ 

"^^  J       r"  4s.  r^-^ 
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Notice  that  since  R,  M  and  Z  behave  like  the  square  root 

of  Y  better  accuracy  is  obtained  by  expressing  their  Y 

2 

derivatives  in  the  form  A  R   /2R.  rather  than  A  R^ , 

<J  J 

Having  found  R^   we  can  now  find  m"?   and  B^"^  from 
J  J        J 

2  2 

5-A-(R^  A+B^)   =   ^  S-[(A-R^  )b'?] 


(5.4) 


2  M^  2 

R^bVIrX)    =   I  [M^A^Or"  +  (60a^)aOr^-  -i(60R^)(A0R'?  )], 


Both  equations  in  (3.4)  represent  linear  differential  equations 

in  the  magnetic  fields.   As  difference  expressions  they 

become  linear  inhomogeneous  equations  for  M.    and  B^ 

J        J 

with  the  inhomogeneity  appearing  as  linear  terms  in 

IVT^  B.,    MV"   and  B.~  .   Thus^  to  an  arbitrary  solution 
J    J    J        J 

of  the  full  equation  we  add  that  solution  of  the  homogeneous 
equation  which  gives  the  desired  value  at  the  boundary. 

Since  both  equations  in  (3.4)  are  of  second  order 
in  Y,    before  we  can  apply  them  we  must  know  the  values  of 
the  fields  at  the  first  two  mesh  points.   The  condition  of 
regularity  allows  us  to  expand  both  in  power  series  about 
the  origin.   Both  equations  also  admit  logarithmic  solutions. 
To  find  these  series  expansions  let  D  =  RM  and  write 

R^(Y,T)  =  2A(A^(T)Y  +  A2(t)Y^  +  A^(T)Y^  +  A^(T)Y   +  ...) 

D  (  Y,T)  =  d^(T)Y  +  d.2(T)Y^  +  d^(T)Y^  +  d^^(T)Y   +  .  .  .  ) 

B  (  Y,T)  =  bQ(T)  +  b^(T)Y  +  b2(T)Y^  +  b^(T)Y^  +  b|^(T)Y^  +  .  .  . 

f(  Y,T)  =  bQ(0)  +b^(0)Y+  b2(0)Y^   ...  =  f^+  f^Y  +  f2Y^+  f^Y^+.  , 
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The  field  equations  are  Integrated  to  give 


(^  ^Y-'y  ~  2"  ^^  ^Y^  ^ 

(3.5)  R  Dyy  "  I  (R^)y° 

(Fy)^  +  I  (F^)y-  >^n2. 


2A(YFy)y  -  A^F 


-^ 


m  R 


0 


When   the    series    are    substituted    into   the    equations    of    (3.5) 
and    the    coefficients    of   like   powers    of   Y   equated    the 
following   set   of   relations    is   found. 

'f^    -    Af^f^      =       A  :n 


(5.6)  { 


4f  f 
1^2 


2  A  f^f^   +    A  f^ 


6f^f^  +  4f2 


8f^f^    +12    f^f^ 


A(5fQf^  +  5   f^f^) 
=      A(f^f,,    +  4   f.f^  +  2   f^) 


0^4 


13 


r 


(3.7) 


4d2    -    Al^ 


m    A  J2l\ 


12i\^d      +  hk^d^    -    A(A^d2    +   2A2d^) 


m  AgA  y2"h 


24A^d^    +   12A2d^    +  4A-,d2    -    A(A^d^    +   2A2d2    +   3A^d^) 


m  A^A  J2l\ 


y 


40  A^d      +   24   A^d^   +  12  Ad      +  4   A^^d^   -    ^(A^d^^    +  2A2d 
+  3  A^dg   +  4  A^^d^)      =      -  m  A^    A  /Tk 
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r 


2A-j_t)-j^  -  M-j^bQ 


2^1  -  ^0 


(5-8) ^ 


8A-j_b2  +  ^Apb^  -  ^(A^b^  +  2A2bQ)   =   Sf^  -  Af^ 
l8A^b^  +  l.'Zk^^   +  6A^b^  -  A(A^b2  +  SA^b^  +  ^A-^b^) 


l8f^  -  Af^ 


32A^b2^  +  24A2b^  +  leA^b^  +  SA^^b^  -  A(A^b   +  2A2b2 
+  3Ajb^  +  4A^bQ)  =  32  f^  -  Af^ 


(3,7)  and  (3.8)  are  used  also  for  the  homogeneous 
difference  equations  by  setting  the  right  hand  sides  equal 

to  zero.   If  we  knew  the  values  of  the  A. 's  we  could  find 

n+1 
these  coefficients.   These  are  determined  from  R.    by 

J 
means  of  the  formula 

2  ^         ^ 


1  T3n+1 
27;  ^j 


4 
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=  A-^jh  +  A2(jh)^  +  A^(jh)^  +  A^(jh)^  +  O(h^) 


in  which  we  let  j  equal  1,  2^,  3  and  4.   Thus 


^  (24r5;+^   -  36  R^^^  +  24r^+^  -6R^+1  )+0(h5) 


(3.9) 


A^h  = 

Agh^  =  ^^  (-26R5+I  +  57R2"*'^ 


42R^ 


A^h^  = 
3 


1     { ^R^+l 


8r 


,n+l' 


+  TR^"^-^  -  2r{J+1  )  +  0(h5) 


5+1   +  11R{J+M-K^(h5) 
2 

2 


,n+l 


There  remains  now  only  Z.    to  find  and  this  is 

J 


obtained  from 
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2  2     2 

(5.10)   5OR-  -^  =  I  5-(|  (  ^  +  ^  )aO  _J_^_^-) 

+  A+(Z^  ,60r^  +  Z^S^R^  J   _ 
^  J-1    J     J    J-1^ 

n+1            64R5+12   2 
with  Z-,     ffiven  by —p-   h  which  is  a  consistent 

representation  for  l/(YRy)  at  Y  =  h. 

n  is  replaced  by  n-1  and  the  process  is  repeated. 

The  difference  scheme  is  accurate  to  second  order 
in  h  and  k. 


b .   The  stability  of  the  system  without  electron  pressure 

In  finding  a  stability  criterion  for  the  system  (3-l) 
we  shall  need  the  definitions 

(A,B)_g  =  h  \__  k.-Q.    ,         WhW^  =    (A,A)_g 

and  the  identities 


,+ 


(A,A-B)^   =   -  (A"A,B)^  +AjBj_^  -  A^_^B^_^ 


(A,A+B)^   =   -  (A~A,B)^  +  ^J-A  -  ^Pl, 

(A^A^B)^   =   -(^°^^B)^+|(Aj_^Bj+AjBj_^-  A^B^_^-  A^_^B^). 
The  energy  of  the  system  (3.1)5  E^  at  time  n  is  defined  to  be 
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E-^+1   =   I      ||5-R^||2    +  ^(IIrH^+b^  ||2+    ||r^+1a+B^+1||2) 


+^(IIa+(rV)||2  +  ||a+(r^+1m^+1) 


+    ||A-(R^M^)||2    +    ||A-(R^+V+1)||2) 


(3.11)  +  ^  (b^b^+\a-(r^%  r''''^^)i 

12"  R^  +  -£+r  ^  ^  (R   +  R     )  Jo- 
in a  manner  similar  to  that  used  to  derive  the 
energy  Integral  (2.54)  we  find  from  (3.1)  and  (3. 11) 
n+1  n  2  2  2 

---K-~  -  h^  "S-i^S  'X-i  -  IT.  («X-i^°«j«"-i"S  ^°«j-i) 

(5.12)  -  ij  By(4:^L-^-_^).   -i^(KX+f5.iMS.i)6V(R5,^) 

+  2^  a  (Zj_j6  Rj  +  Zj5  Rj_-]_)  =  0(1^  )  • 

A  derivation  of  (3. 12)  is  given  in  Appendix  2.   This  equation 
says  that  the  total  change  in  the  energy  of  the  plasma  is 

equal  to  the  energy  gained  at  the  boundary  plus  terms  of  the 

2 
order  k  ,    which  is  the  accuracy  of  the  difference  method. 

n 
If  E   is  a  valid  energy^  that  is  if  it  is  positive,  the 

difference  scheme  will  "be  stable. 

A  sufficient  condition  that  E^  be  positive  is  given  by 
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(5.5)   ,  (^  „n^-i  H-  ^4  (| .  |^))A°(Rf  +  Rr''" 

j   J 

+  Z^    Z""^        2        2 

The  last  term  in  (3.5)  is  indeterminate  as  written  and 
is  found  from  the  series  expansion  of  R. 
Let 

L^  -   ^  ^,^r\  ^;S   .  !|;J)  and  I-  .  ^  B>-^  llL^.L^^), 


Since 


-.•_n     J        J        J  1=1     '^    J     <J   -^  <J        J 

-|LgA-(R5;\R-^')+|L5.,A-(Rf-R5-^') 


the  condition  (3.5)  can  be  written 

l"^        ,  ,  -  -  ^  ^  ...   ... 


tz    [(5-R")'  +  I-A-(Rf  +  R^-l')]  -  I  L2A-(Rf  +  R^-l") 


The  last  three  terms  arise  from  the  presence  of  the 
boundaries.   Their  sum  is  positive  and  does  not  affect  the 
main  part  of  the  stability  criterion^  the  summation.   We 
can  write  this  summation  in  the  form 
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(3.15)  ZZ  (|  ^"^5)^+  |(5"R-_i)^+  ^^"(l^^+R^"^^] 

The  last  term  is  positive  and  5~r!^  is  zero.   Writing 

2     ^_^2  ^_^   _  ^    ^ 

r'^   as  R^    +  (R^  +  R.~  )5  R^k  we  see  that  each  term 

J  J  J         J  J 

in  the  series  in  (3. 15)  will  be  positive  if 

(3.16)  I    (5-R^)2+  1(5-r'?   )2+  r^A-(2R^-l%(R^+R^-l)5-R^k)>0 
for  0  <  j  <  J.  A  sufficient  condition  that  (3.I6)  holds  is 

(3.17)  i(5"R?)^+ i(6"R'?_i)^+  2A-r']-'^  >   |i^a-((r^+r^-Vr'?)|i^. 

^J^J-*"  JJ  J         JJ         J 

However^  the  right  hand  side  of;  this  satisfies 

klrJ'A-((R;j+R|J-i)6-R'}l  iJi|r^|{|{R5+R5-^)6-R;}| 

^l(H^.-H-i)6-Rj.J)<|^rf({Rn.Rn-l^ 


-.n    ,  „n-l\2N  ,  1  /5---p,nN2  ,  1  /o-^n  \2 


Upon  substituting  this  into  (3.I7)  we  arrive,  after  cancel- 
lation,  at  the  stability  criterion 


■'  It  a-r"-^' 

J 

Since  the  sound  and  Alfven  speeds  are  given  by 
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C?  =  -^  and  C?  =  MB!+M!i 

the   stability   cirterion    (3-l8)    specifies,    discoujiting  the 

centering  of   the   quantities    in   the   mesh,    that 

,hs2        12  12 

^-^)      "^  If  ^A    +  12   ^S- 

A  stability  check  is  performed  at  each  step  in  the 
numerical  calculations  using  (5.I8).  Also,  the  energy 
equation  {2.'jh)    is  integrated  numerically  as  a  check  on 
the  accuracy  of  the  calculations.   The  two  parts  of  the 
integral  still  agree  to  better  than  five  parts  in  a  thousand 
at  time  T  =  20. 

We  are  now  in  a  position  to  justify  the  use  of  the 
system  (3.1)  in  solving  this  problem  rather  than  using  (2.48) 
as  it  stands.   Whereas  the  equation  for  B  lends  itself 
readily  to  incorporation  into  the  discrete  energy  equation, 
it  was  found  that  the  equation  for  M  did  not.  The  equation 
for  M  that  appears  in  (3.1)  was  rigged  in  such  a  manner  as 
to  make  all  troublesome  terms  cancel  in  the  final  addition. 
For  this  convenience  we  pay  a  small  price  in  the  complexity 
of  the  equation.   The  trouble  with  leaving  Z  as  l/YRy  appears 
when  the  energy  equation  is  derived  in  discrete  form  to  prove 
stability.  The  object  is  to  make  the  change  in  energy  in  one 
time  step  proportional  to  k^.   However,  we  find  that  in  the 


derivation  we  use  differential  identities  in  which  Y  appears. 

o 
Such  terms  can  be  correct  only  to  0(h  ).   Instead  of  having 

the  change  in  energy  0(k  )  we  could  obtain  it  only  to  0(h  +k'^) 
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and  thus  not  sufficient  for  proving  stability.   By  using 
the  differential  equation  to  define  lARy  ^^  have  eliminated 
this  explicit  dependence  on  Y  and  have  forced  the  energy 
equation  into  the  desired  form.   Of  course,  we  have  had 
to  introduce  another  variable  and  another  equation. 

c .   Numerical  simulation  of  the  pinch  with  pressure 

The  inclusion  of  the  electron  pressure  causes 

considerable  complication  of  the  numerical  method  used  to 

solve  the  problem.  This  arises  primarily  from  the  time 

derivatives  which  appear  on  the  right  hand  sides  of  the 

field  equations.   While  it  seems  possible  to  solve  the 

new  system  numerically  we  have  not  been  able  to  devise  a 

_2 
stable  method.   However,  by  considering  terms  of  0(a   ) 

negligible  a  stable  and  greatly  simplified  method  can  be 

used. 

Again,  scaling  the  Y  variable  by  the  parameter  A, 

(2.48)  becomes 

J.  (9L-h—\   _  ^^  Y  ,  /g  P  \   _   a   A 

2YRy5^Y"    r3-  -^  ^  ^  ^Y  -j,3i^2p232  = 

(3.19)       ^     oTf  §h2tA''\      • 

%  "12   ^  Hg   ^YY  '   ^T  ~  ~  ^^"^     H^     Yy' 

Hg   +   t;  RyBg  -  -^  (RB2)Yy  '>    ^^    =   T   ^Y^^"  ~2^^  ^^^Y^ 
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The  equations  for  Hp  and  H^  in  (5.I9)  can  be  written 
in  the  form 

R  7v     -  S         -  S 

"^T  2  ^      Hp   T  ^^ 

(3.20)  2     ^  2 

^3^  =  -  27^  La  ^— — ^T  -  "  ^  ^3/  °^  "~J~  ^2^JyY' 

However,,  since  Hp^  and  H^m  are  0(a~)  these  equations  give 
.2    -r.2  2 


12~  ^^"^TYY  +  0(°^   ) 


^2 

^^    --     -2 

(3.21)  2 

/-/,  T,2    _  H^S  p 

%  =  "  277f"  ^^"   h2"^TYY  +  0^°^"  )  • 

Neglecting  the  higher  order  terms  these  equations  can  be 
integrated  with  respect  to  T  using  the  initial  conditions 
given  in  (2.53).  Further  simplification  can  be  achieved  if 
we  notice  that  H^,  and  H^^  as  well  as  all  Y  derivatives^ 
remain  within  0(a~)  of  their  initial  values^  denoted  by 
Hppj  and  H^p^.   Since  Hpj  H^  and  S  are  initially  constant 
(3.21)  when  integrated  yields 

2  2 

\  =   H20  +  g-  g^  (a-S^)^^  =  H2  +  0(a-  ) 

(3.22) 

u    u     64a2  H  q   _  2      ^     ^  CM    -^\ 

"3  =  ^30  -  27Tf-  Zt  ^"^   2  )yY  =  H^  +  0(a   )  . 

^20 

Since  we  are  solving  this  problem  to  first  order 

in  a~  we  see  that  it  is  only  necessary  to  find  P  and  S  to 
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this  order  since  whenever  they  appear  in  (3.19)  they  are 

multiplied  by  a  factor  a.~ .      We  therefore  introduce  an  error 

2 
of  0(a~  )  by  using 


2  2 

Ho  ^20   '        RRy^    ^30      ^2 


(-.    o^)    ^6  _  1  ^r   T   S5Y  ^2    .   -6  _   1  h~   ^  ^3Y  .^20^2 


in  place  of  P  and  S  in  (5.I9)  and  (5.22).   We  thus  obtain 

from  (3.19)  the  following  system  which  differs  from  it  only 

_2 
by  terms  of  0(a~  ). 

1         2    1  1 

^  .f    q+ ^    2?v:x"^Y  _^  /a"P^^        a'A^ 

+  A( ^) ^-  +  (-^— ) ^^;^^^  =  0 

2YRy'^  Y    R^      ^y   1    R(RRyPS)'^ 

^2   "   ^20  +  T2IU7  ^°^  ^^YY 
(3.24) 

H    -   H      ^^   --50  (a--.) 
3  ~    50  ~  27  7r  „2      ^      ^YY 

^20 

(RB2)Yy  -  ^^2    +  ^^^2   ^   ° 

(R^B   )  -  ARRyB   +  A^H^  =   0  . 

By  setting  B^  and  H^  equal  to  zero  in  (3.24)  we 
obtain  the  equations  used  in  the  numerical  simulation  of 
the  z  pinch.   For  the  theta  pinch  the  system  must  be 
modified  to  correspond  to  (2.5O).  When  this  is  done  we  obtain 
a  new  equation  for  H^  given  by 


fV* 


25 


^3   ~   ^30  ~  tSTT   H^„  ^r,2^2  ^YY 

3U   r  Ky 
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Bp  and  Hp  are  set  equal  to  zero  in  the  remaining  equations. 

The  numerical  solution  of  this  system  is  now  not 
unlike  that  of  the  one  in  which  the  electron  pressure  is 
absent.   Assuming  we  know  all  the  variables  at  two  time 
steps  we  find  R  at  the  next  time  level  from  the  first 
equation  in  (3.24).   We  then  find  B^  and  B  using  the 
initial  values  of  H^  and  H^.    This  will,  of  course,  be 
correct  only  to  0(a")T5  but  that  is  sufficient  to  find  the 
new  values  of  P  and  S  since  they  are  desired  only  to  this 
accuracy.   With  the  new  P  and  S  we  find  H^  and  H^  and  from 
these  we  obtain  corrected  values  of  B„  and  B^ .   The 
procedure  is  now  repeated  to  advance  in  time. 

Although  we  have  no  proof  of  the  stability  of  this 
method  we  may  infer  from  the  well  behaved  results  and  the 
results  of  similar  problems  for  which  stability  has  been 
proven  that  the  method  is  stable. 
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CHAPTER  IV.   NUMERICAL  SIMULATION  OF  THE  PINCH 

The  systems  (2.48)  and  (2.50)  with  the  initial  and 
boundary  conditions  given  in  (2.55)  and  (2.55a)  are  solved 
numerically  by  the  methods  described  in  Chapter  III.   The 
results  of  those  computations  are  presented  in  this  chapter. 

a .   The  pure  theta  pinch 

System  (2.50)  representing  a  theta  pinch  is  solved 
under  the  following  conditions.   To  a  column  of  plasma  of 
radius  ten  satisfying  the  equilibrium  conditions  given  in 
(2.53)  we  raise  B^  linearly  at  the  boundary  to  double  its 
initial  value  in  a  time  Tp.   Thereafter  B^  is  held  at  two. 
The  parameter  To  regulates  the 'strength  of  the  compression 
by  specifying  the  length  of  time  the  plasma  feels  the  full 
strength  of  the  field.   Figures  1,  2  and  5  show  the  density 
profile  of  the  plasma  column  for  pressure-free  cases  in 
which  Tp  equals  10,  7  and  2  respectively.   In  all  three 
we  see  a  train  of  waves  being  generated  by  the  boundary  and 
travelling  inward  ahead  of  the  incoming  plasma.   Upon  reach- 
ing the  origin  the  density  rises  to  several  times  its 
initial  value.   The  waves  are  then  reflected  outward  where 
they  interact  with  the  remaining  tail  of  the  wave  train. 
In  Figure  1,  showing  the  mildest  compression,  the  first  wave 
reaches  the  origin  at  approximately  T  -   10.1  and  the  density 
rises  to  a  value  of  about  5.6.   Since  our  scaling  velocity 
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is  the  Alfven  velocity  and  our  initial  radius  is  ten  we 
see  that  the  wave  travels  at  approximately  the  Alfven  speed, 
as  a  linearized  theory  would  predict.   In  Figure  2.,    for 
which  T^  =  7,  the  first  wave  reaches  the  origin  at  T  =  9.8 
and  the  density  rises  there  to  9.'^.    In  the  much  stronger 
compression  produced  by  T^  =  2  and  shown  in  Figure  5  the 
disturbance  arrives  at  the  center  at  T  =  7-2  causing  a 
density  rise  in  excess  of  JO.   Notice  that  in  Figure  3  when 
the  first  wave  is  reflected  from  the  origin  it  draws  with 
it  most  of  the  gas  so  the  density  falls  rather  sharply  at 
the  boundary.   When  the  second  wave  reaches  the  center  at 
T  ^  11  the  density  is  seen  to  again  rise  quite  strongly. 
Finally,  we  see  the  first  wave  has  traveled  back  through 
the  plasma  and  has  been  reflected  at  the  outer  boundary. 
Since  the  plasma  as  a  whole  is  now  moving  outward  the 
reflected  wave  gains  in  strength.   At  T  =  l4  a  good  part 
of  the  plasma  is  concentrated  in  this  wave. 

Figure  4  shows  the  effect  of  pressure  on  the  motion. 
Again  taking  T^  =  7  but  now  setting  a"^  =  a"  =  .04  we  see  by 
comparing  Figure  4  with  its  pressure-free  counterpart  Figure  2 
that  the  pressure  causes  two  obvious  changes  in  the  motion. 
In  addition  to  the  initial  Alfven  speed  of  unity  we  now  have 
an  initial  sound  speed  which  equals  J^a     +  Ja"  ,    or  /TW  . 
Thus  instead  of  the  disturbance  taking  nearly  10  time  units 
to  reach  the  center  it  now  takes  about  R^/  sJz.-vQ,^   =  10/  >/L+. 2^ 
%  8.96.   In  the  actual  numerical  simulation  the  first  wave 
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reaches  the  origin  at  about  8.82.   In  addition  to  the 
increase  in  the  speed  of  propagation  we  observe  another 
effect,  the  tendency  of  the  waves  to  be  more  sharply  peaked. 
The  pulses  are  no  longer  rounded  as  they  were  in  Figure  2 
but  have  a  steeper  front  and  sharper  peak.  Wnen   the  pulse 
reaches  the  origin  the  density  rises  above  the  value  IJ , 
nearly  twice  that  of  the  pressure-free  case  with  its 
broader,  more  rounded  wave.   The  effect  of  increasing  a 
and  a"  is  to  cause  the  wave  fronts  to  steepen  until  the 
density  becomes  discontinuous. 

The  axial  magnetic  field  is  approximately  equal  to 
the  density,  differing  from  it  only  in  that  the  pulses  tend 
to  be  weaker  and  lag  slightly  behind  those  of  the  density. 
These  deviations  from  the  one-fluid  fields  are  produced  by 
the  terms  in  H^  involving  derivatives  of  B^. 

Kilb  (12)  has  obtained  similar  results  in  theta 
pinch  calculations. 

b .   The  pure  z  pinch 

Corresponding  to  equations  (2. 51)  and  (2.52)  is  a 
compression  in  which  the  only  field  present  is  Bp.   Unless 
we  take  a  non-constant  initial  pressure  the  only  equilibrium 
is  one  in  which  Bp  is  identically  zero.   If  we  start  the 
problem  with  no  magnetic  field  we  have  to  define  a  new 
scaling  parameter  in  (2.4).   We  take  for  this  the  final 
value  of  the  magnetic  field  at  the  boundary.   Thus,  the 
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boundary  conditions  on  Bp  for  this  problem  are 


(4.1)  Bp(T,Yj  .  ^^f  ^ 


^1    .   T^  <  T  . 

¥e  find  that  with  this  pinch  it  is  impossible,  even 
with  very  mild  compressions ,  to  obtain  solutions  which  do 
not  break  almost  im.mediately .   The  reason  for  this  is 
evident.   Without  a  magnetic  field  in  the  plasma  the  only 
mode  of  propagation  initially  present  is  the  pressure. 
Thus,  until  the  field  penetrates  into  the  plasma,  the 
behavior  is  similar  to  an  ordinary  fluid,  that  is,  disconti- 
nuities tend  to  develop  upon  compression.   Since  we  are 
concerned  in  this  paper  with  well  behaved  flows  we  shall 
not  pursue  the  investigation  of  the  pure  z  pinch  but 
rather  examine  a  hybrid  pinch  in  which  a  magnetic  field  is 
initially  present. 

c .   A  hybrid  z-theta  pinch 

In  this  pinch  we  start  with  the  same  equilibrium  as 
in  the  theta  pinch,  both  N  and  B^  equal  to  one.   In  this 
case,  however,  the  motion  is  produced  by  raising  Bp  at  the 
boundary  while  imposing  the  condition  that  B-^y  vanishes 
there.   With  this  condition  on  B^y  no  axial  component  of 
the  field  enters  or  leaves  the  plasma  at  the  boundary  and 
we  shall  show  that  the  total  flux  of  B^  in  the  gas  is  constant 
to  0(a~).  The  time  derivative  of  this  flux  is 
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Rm(t) 


■M 


■M 


_d_ 

dT 


0 


rB^(t,r)  dr  =  1^  ^ 


0 


0 


+  R-B 


1-  r 

0 


r  RRyB^  dY  =  J^-  1  H^  +  (R^B^y)^^! 


Y 


M 


3Y 


1   . 


We  have  used  the  definition  of  H^  given  in  (2.47).   By 

virtue  of  the  boundary  conditions  the  second  term  vanishes 

while  the  integral  can  be  evaluated  using  (2.48)  to  yield 

Y, 


•M 


■M 


0 


H^   dY  = 


64a~ 
27  f' 


64a_ 
27  TT 


( 


2 


0 


H 


^  )yy'^y 


2 

H  S 
,3   T, 


Y=Y 


H 


2  ^Y 


M 


Y=0 


Thus  J 


(4.2) 


Rm(t) 


'M 


dT 


64 


a 


J  rB^(t,r)  dr  =   -    ^^ ^ 
0 


2 
(-\— )^ 


Y=Y 


H 


2   ^Y 


M 


Y=0 


giving  us  the  desired  result 


B 


P  is  raised  linearly  at  the  boundary  from  its 
initial  value  of  0  to  a  value  of  .6  in  a  time  T„  =  10. 
Thereafter  it  is  held  fixed  at  .6.   To  allow  the  computa- 
tions to  proceed  without  the  solution  breaking  we  take 
the  initial  pressure  to  be  quite  small:   a  =  a  =  .001. 
Because  of  the  high  densities  encountered  in  this  case  the 
pressure,  even  with  this  small  initial  value,  becomes 
significant. 
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In  Figure  5a  we  observe  that  this  pinch  takes  longer 
to  develop  than  a  theta  pinch.   In  the  latter  the  disturbance 
has  reached  the  origin  at  T  =  10  while  in  this  pinch  it  has 
propagated  only  halfway  to  the  origin.   This  delay  occurs 
because  the  magnetic  field  enters  the  momentum  equation 
quadratically  and  whereas  the  compressing  field  was  initially 
one  in  the  theta  pinch  it  is  initially  zero  in  this  pinch. 
Thus,  it  takes  a  while  before  sufficient  Bp  is  present. 
Once  the  motion  does  begin,  however,  it  is  much  stronger 
than  anything  observed  in  the  theta  pinch.   One  reason  for 
the  rapid  compression  is  that  as  the  density  rises  near  the 
bo^ondary  so  does  B^  and  this  aids  in  the  initial  compression 
stages.   A  more  fundamental  difference  is  that  in  the  z  pinch 
there  is  no  equilibrium  state  nearby  about  which  the  plasma 
column  can  oscillate.   As  we  saw  in  the  theta  pinch,  the 
radius  does  not  continue  to  diminish  but  rather  oscillates 
about  a  near-equilibrium  state.   Indeed,  with  any  dissipa- 
tion present  the  plasma  would  eventually  approach  an  equili- 
brium of  constant  density  and  constant  magnetic  field.   In 
the  z  pinch,  though,  the  plasma  is  not  close  to  an  equilibrium 
state  and  so  the  radius  continues  to  fall.   To  stop  the 
collapse  it  would  be  necessary  to  take  into  account  the  jump 
conditions  at  the  plasma-vacuum  boundary  where  the  pressure 
will  eventually  become  strong  enough  to  balance  the  field. 
Since  we  have  ignored  this  effect  our  solution  would  soon 
become  singular.   These  considerations,  however,  do  not  affect 
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the  validity  of  the  problem  in  the  range  illustrated  in 
Figure  5. 

Referring  now  to  Figure  5t)j  we  see  that  by  T  =  12 
the  radius  of  the  column  has  shrunk  to  ^-'..9.      At  T  =  12.8 
the  first  pulse  reaches  the  origin  where  it  peaks  to  a 
value  greater  than  35.   At  T  =  I3  this  pulse  has  been  reflected 
outward  toward  the  boundary  which  is  now  at  a  distance  of  only 
5.5  away  from  the  origin  and  still  moving  inward.   Finally, 
by  T  =  l4  the  radius  has  fallen  below  2  meaning  that  the 
average  density  is  greater  than  JO.   During  the  final  stages 
of  the  collapse,  local  density  values  of  near  100  are  observed. 

The  noteworthy  characteristic  of  this  flow  is  that 
very  strong  compressions  and  lq,rge  density  rises  are  obtained 
without  initiating  the  formation  of  discontinuities  even 
though  the  value  of  B^  which  must  be  imposed  on  the  boundary 
to  satisfy  the  condition  B^y  =  0  reaches  above  20.   In  a 
theta  pinch  it  would  be  impossible  to  increase  B^  to  this 
value  so  quickly  without  the  solution  breaking  down.   Part 
of  the  reason  for  the  well  behaved  nature  of  the  solution 
is  that  the  value  of  B  at  the  boundary  is  adjusted  in  such 
a  way  that  there  is  no  force  produced  by  this  component  at 
the  boundary. 

If  the  initial  pressure  is  taken  to  be  larger, 
discontinuities  in  the  density  develop.   Even  with  fairly 
substantial  initial  pressure,  though,  the  solution  remains 
well  behaved  until  large  densities  are  reached.   For  instance, 
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with  an  initial  pressure  of  a  =  a~  =  . 05  the  radius  still 
falls  below  4.5  without  the  solution  breaking.   If  an 
artificial  viscosity  is  added  it  appears  that  the  smoothness 
of  the  solution  can  be  maintained  and  the  shock  resolved. 
It  should  be  mentioned  that  the  ratio  3  =  (P-^-j^  +  ^11^/2"  ^ 

remains  small  even  though  the  pressure  increases  to  many 

2 
times  its  initial  value  because  B^,  like  the  pressure^  xs 

approximately  proportional  to  the  density  squared. 

d .   Other  flows 

When  two  components  of  the  magnetic  field  are  present 
in  the  plasma  we  observe^  by  special  choice  of  the  boundary 
conditions^  unusual  phenomena  which  arise  from  the  distor- 
tion of  one  component  by  the  compression  initiated  by  the 
other  component.   These  effects  depend  strongly  on  the 
initial  geometry  of  the  fields  and  can  become  so  strong  as 
to  prevent  the  numerical  computations  from  proceeding. 

We  take  the  pressure  to  be  zero  in  all  these  cases. 
The  equilibrium  is  determined  by  the  condition  that  forces 
produced  by  the  fields  balance.   Thus^  the  appropriate 
initial  conditions  are  given  in  (2.55a).   From  this 
equilibrium  one  of  the  magnetic  fields  is  raised  at  the 
boundary  while  the  value  of  the  other  field  at  the  boundary 
is  held  fixed.   It  is  this  condition  which  produces  the 
unusual  flows. 

In  Figure  6  we  show  the  density  profile  which  results 
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from  raising  Bp  "while  keeping  B^  equal  to  zero  on  the 
boundary.   The  initial  fields  are  shown  in  the  upper  left 
hand  corner  of  the  graph.   Figure  6  should  be  compared  with 
Figure  1  in  which  B^  was  raised  at  the  same  rate.   Two 
differences  are  seen  immediately.   First;,  the  waves  are 
weaker  and  broader.   The  most  striking  difference,  though, 
is  the  sharp  rise  in  density  at  the  boundary  which  occurs 
in  Figure  6  but  is  absent  in  Figure  1.   This  effect  is  more 
pronounced  in  Figure  7  in  which  more  B  and  less  B„   is 
present  initially.   Again,  the  total  field  strength  at  the 
boundary  starts  at  one  and  is  raised  linearly  to  two  and 
then  held  fixed.   If  we  start  with  the  equilibrium  in  which 
no  Bp  is  present,  that  is  as  we  did  in  the  pure  theta  pinch, 
but  this  time  raise  Bp  we  find  that  for  the  same  strength 
compression  the  density  at  the  boundary  becomes  so  large 
so  fast  that  the  numerical  computations  fail. 

To  understand  this  phenomenon  we  must  look  at  the 
geometry  of  the  fields  and  the  imposed  boundary  conditions. 
In  Figure  11  we  have  shown  the  density  and  two  magnetic 
fields  for  the  problem  considered  in  Figure  Y,    that  is, 
starting  from  the  equilj.brium  in  which  Bp  and  B^  are  equal 
on  the  boundary  we  raise  Bp.   As  Bp  is  increased  at  the  edge 
the  plasma  moves  inward  and  the  density  near  the  boundary 
rises.   B^,  behaving  like  the  density  also  rises  near  the 
boundary.   The  value  of  B^  at  the  boundary  is,  however,  held 
at  its  initial  value.   This  creates  an  increased  B^  inside 
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the  plasma  which  tends  to  inhibit  the  inward  motion  of 
the  gas  at  the  boundary.   Thus  the  gas  near  the  boundary 
is  caught  between  the  rising  Bp  field  at  the  boundary  and 
the  rising  B-,  field  inside.   This  gives  rise  to  the  strong 
increase  in  density  at  the  boundary.   These  events  can  be 
seen  in  Figure  11.   If  the  amount  of  B^  present  is  not  too 
large  the  compression  wave  can  overcome  the  barrier  and 
propagate  into  the  plasma.   When  this  happens  the  B^  peak 
also  propagates  inward,  temporarily  relieving  the  pressure. 
This  can  be  seen  in  Figure  6  at  times  5  and  9  and  in 
Figure  7  at  times  6  and  9-   In  the  case  where  there  is  no 
Bp  present  initially,  though,  the  density  has  already  risen 
past  computable  values  before  the  B^  barrier  can  propagate 
inward  and  relieve  the  pressure.   Because  any  compression 
waves  that  do  propagate  inward  must  first  overcome  this 
resistance  they  are  weaker  than  their  counterparts  in  the 
pure  theta  pinch. 

There  are  other  factors  which  contribute  to  the 
rapid  rise  in  density  produced  when  Bp  is  used  to  compress 
the  plasma.   These  are  related  to  the  means  of  propagation 
which  are  available  to  a  disturbance.   Since  there  is  no 
pressure  in  this  problem  the  only  means  available  is  through 
magnetic  compression  waves.   These  are  produced  by  the 
magnetic  fields  acting  as  pressures  in  the  momentum  equation 


:^.3)      RrpT  +  I  RlBp  +  B^)y  =    R(B^y)^  -  V2 
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The  field  terms  on  the  left  are  the  magnetic  pressure  terms 
while  those  on  the  right  are  related  to  the  cylindrical 
geometry  and  are  in  general  much  smaller  away  from  the 
origin.   If  we  attempt  to  compress  the  plasma  using  Bp 
and  start  from  an  equilibrium  in  which  no  Bp  is  present 
we  see  that  since  B^  appears  in  (4.5)  only  quadratically 
it  is  very  inefficient  as  a  means  of  wave  propagation. 
Since  B-,  is  held  fixed  at  the  boundary  it  cannot  help  much 
in  propagating  waves  inward  from  the  edge.   Once  inside 
the  plasma,  though,  B^  becomes  an  effective  means  of  wave 
propagation  since  it  is  no  longer  restricted.   In  fact, 
near  the  origin  it  becomes  the  only  means  of  wave  propaga- 
tion since  Bp  must  vanish  there.   If  we  provide  another 
means  of  propagation  by  adding  ion  pressure  to  the  problem 
the  peaking  is  reduced. 

There  is  another  boundary  phenomenon  which  occurs 
when  both  fields  are  present.   If  we  again  start  from  that 
equilibrium  in  which  Bp  is  one  at  the  boundary  while  B^ 
vanishes  there  and  raise  B^  instead  of  Bp  we  obtain  a  very 
different  result.   Whereas  when  Bp  was  raised,  as  illustrated 
in  Figure  5,  the  density  rose  at  the  boundary,  when  B^  is 
raised  it  falls  1   This  is  shown  in  Figure  8.   Notice  that 
this  effect  is  quite  strong  with  the  density  at  the  boundary 
falling  to  about  one-tenth  of  its  initial  value. 

The  explanation  for  this  behavior  again  lies  in  the 
very  special  geometry  of  the  initial  fields.   Initially, 
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to  balance  the  force  produced  by  the  positive  constant 
slope  of  Bpj  B^  has  a  negative  slope.   B^,  however, 
contributes  two  terms  to  the  equilibrium^  the  magnetic 
pressure  B^B^y  and  the  centrifugal  term  (B^y)  .   In 
general^  the  pressure  terai  is  much  larger  but  in  this  case 
we  have  made  B-,  vanish  at  the  boundary  and  so  the  centri- 
fugal  term  alone  must  provide  the  force  needed  to  maintain 
equilibrium.   This  gives  rise  to  a  large  negative  slope  for 
B^  at  the  boundary.   When  B^  is  raised  at  the  edge  two 
things  happen.   Where  before  there  was  no  B^j  and  thus  no 
magnetic  pressure^  there  now  is  magnetic  pressure  and  since 
the  slope  of  B^  is  large  and  negative  the  effect  of  the 
pressure  is  to  drive  the  plasma  outward.   Thus^  even  though 
we  are  increasing  B  ^  the  effect  at  the  boundary  is  to  push 
the  gas  away  from  the  center.   There  is  another  effect 
causesd  by  the  increasing  B^  which  also  contributes  to  the 
density  drop.   As  the  column  as  a  whole  contracts  under  the 
increasing  B  field  the  density  inside  rises^  and  so^ 
therefore,  does  Bp.   Since  Bp  is  fixed  at  one  at  the  boundary 
the  rise  of  Bp  inside  the  plasma  has  the  effect  of  diminishing 
the  slope  of  Bp  at  the  boundary.   Thus,  the  ability  of  B^  to 
contain  the  plasma  is  diminished  at  the  edge  as  the  B^  field 
grows.   The  density  at  the  boundary  thus  falls.   These 
developments  are  illustrated  in  Figure  10.   Notice  in 
particular  how  the  slope  of  Bp  is  diminished  from  its  positive 
initial  value  and  hoW;,  by  T  =  10,  it  has  become  negative,  thus 
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contributing  to  the  outward  accelration  of  the  plasma  at 
the  boundary. 

This  density  drop  is  a  direct  consequence  of  the 
fact  that  B  vanishes  initially  at  the  boundary.   For  when 
B^  is  initially  present  in  any  non-trivial  amount  this 
drop  does  not  occur.   Figure  9  shows  the  density  profile 
produced  by  raising  B^  from  the  initial  state  in  which 
both  fields  are  equal  at  the  boundary.   The  drop  in  density 
has  disappeared  and  we  now  even  see  a  slight  rise  in  the 
density  at  the  boundary  caused  by  the  same  mechanism 
described  earlier  with  the  role  of  the  fields  reversed. 
By  increasing  the  initial  value  of  B,  still  further  we 
approach  the  pure  theta  pinch  of  Figure  1. 

Since  the  boundary  effects  discussed  in  this  section 
are  primarily  local  phenomena  they  depend  on  the  initial 
radius  of  the  column  only  insofar  as  it  affects  the  local 
field  geometry.   Thus,  these  phenomena  are  not  associated 
with  cylindrical  compression  but  arise  whenever  one  attempts 
to  compress  a  plasma  which  contains  a  magnetic  field  which 
is  held  fixed  at  the  boundary. 


-75- 


APPENDIX  1.   Characteristics 

In  this  section  we  derive  the  characteristic 
equations  given  in  (2.54).   We  examine  first  the  simple 
case  in  which  there  is  no  pressure.   (2.48)  then  becomes 

R^^  -  R(B^y)^  +  B2(RB2)y  +  RB^B^Y  =  ^ 
(Al.l)  [(RB^)^  -  R^Bgl^  =  0 

[(R^B^y)y  ~  ^V^^T  =  °  • 

From  the  second  two  equations  we  immediately  obtain  the 

o 
characteristics  (dY)   =  0.   After  integrating  we  can  write 

the  equations  as 


xVprp     +     .  •  •      —     w 
(A1.2)  RB2YY     +     B^Ryy     +      ...      =      0 

2 
R  ^-jYY  +  •  •  •  =  ^ 

where  the  dots  refer  to  terms  with  lower  order  derivatives. 
The  last  equation  gives  the  characteristics  (dl)   =  0. 
¥e  are  thus  left  with 

Rmm  +  . . .  =  0 
(A1.3)  ^^ 

RB2YY  "*"  ^2'^YY  +  •  •  •  =  "^ 

which  can  easily  be  shown  to  contribute  the  characteristics 

p 
(dT  dY)   =  0.   Thus  the  full  characteristic  equation  is 

4   4  2   ^ 

dT  dY  =  0.   For  the  pure  z  pinch  this  becomes  dT  dY  =  0 
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2    ^ 
and  for  the  pure  theta  pi^nch  we  obtain  dT  dY  =  0. 

Tlie  inclusion  of  ion  pressure  changes  the  first 

equation  to 

while  leaving  the  others  unchanged.   After  integrating 
and  obtaining  the  characteristics  (dY)   =  0  the  system 


becomes 


+2 

■f^mm  ^T  YY  •    •    •       —      U 


(A1.4)  RB2YY   +  B^RyY   +    ...    =    0 

2 

R        B-7yy.  +  .      .      .  =  0  . 

p 

From  the  last  equation  we  again  obtain  (dT)   =  0  while 
the  first  two  yield  dT2(dY2  -  c^  dT^)  =  0.   Thus  the 
characteristic  equation  for  the  system  is 

2 

-dT^     "L 


dT^  dY^  [{^f   -  cf  ]  =  0  . 


The  characteristic  equation  for  both  the  z  and  theta 

pinches  in  this  case  is  given  by 

p 
dT^  dY  [(g)2  -  cj  ]  =  0  . 

We  consider  now  the  problem  in  which  both  the 
electron  and  ion  pressures  are  included,  examining  first 
the  pure  theta  pinch.  The  system  (2.5O)  can  be  written 
in  the  abbreviated  form 
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2        -  2 
^TT  "  ^L  ^Y  ~  ^"2  ^YY  "^  "rT  ^  Y"^  '  '  '  "  ° 
-  2 

^     TYY      ?ry~P 

(A1.5)  R[  --^-^^  -      ^^^.^TYYY^  +  •••   =  ^    '    ^"^^""^ 

Y  Ry*~ 

4      ^5  ~  ^^5Y   2 
^Y   =   ^R^HZ    ^  ^5YYY  "^  •  •  • 

From  the  first  two  equations  we  obtain  the  characteristic 

equation 

2     2 
dT^  dY  [dY^  -  (c+  +  c£  )  dT^ ]  . 

However,  we  wish  to  solve  for  B  and  R,  not  P  and  R 

so  the  characteristic  condition  for  the  pure  theta  pinch 

becomes 

dT   dY  [dY   -  (cj  +  c£  )  dT  ]  =  0. 

The  system  (2.52)  representing  the  pure  z  pinch 

can  be  written  in  the  form 

2        -  2 
^TT  "  ^L  ^Y     ~2~   \y  "^  R;;  ^  Y  ■^-  •  •  "  ° 

Ky  1 

(A1.6)  S    rpYY+    •..    =    0     but   now    P^   is    given   by 

4  Ry'^2         1 

^^      "      '    (RR    )2      JJ2    (^^2YY    "^  B2\y''     • 

With  this  definition  of  P  we  see  that  the  third  derivatives 

2 
of  R  now  appear  in  the  first  equation,  making  the  P   term 

dominate.   However,  we  see  that  the  dominating  parts  of  the 

first  two  equations  are  equal.   We  must  therefore  eliminate 

P  in  the  first  equation  using  the  second  to  obtain  an  equation 


-76- 


with  independent  dominating  terms.  This  is  done  in  the 
following  manner.   Using  PS  =  1/(rRy)  the  system  becomes 


+       a~P        a~  /  1   2    i^Y    Y 
^TT  "  ^L  ^YY     r~2"  ^YY  ~  RT  ^W^'^    ^"RT"  "^"72^^-  '■    =   ^ 


Ry 


Y   ^Y     'T    S' 


2         2Y   2      2YY   2 
(AI.7)  S  rpYY"  2  -j^  S  ^Y~  ~lj —  S  rjn  +.  •  .  =0  where 


g 
Sy  =   p  Py  +  . . .   and  Py  is  given  above, 


We  now  differentiate  the  first  equation,,  remembering  that 

this  adds  characteristics  which  must  later  be  eliminated. 

This  gives 

2         -  2 
+         a  P 

\tTY  "  ^L  \yYY     r~2  ^TYYY. 

Ry 


«:  (_1_)2  (^^YY    ^  TYYx         _   0 
Ry  ^RRyS^   ^   Ry        g2  ^  ^  • • • 


When  we  use  the  second  equation  to  eliminate  S^yy  we  are 
left  with  the  system 


(A1.9) 


2 

R        (c+^+c-^)T^  2«-p^   ^2    c^i^   +  =0 

TTTY  "  ^  L      L  ^^TYYY      Ry     H^   TY   Ry   "H^2YY  *" 

2 

T'W  "*   •  •  •   ^   w  • 


These  two  equations  are  independent  and  to  find  the  charac- 
teristics we  use  the  definitions  of  Hp  and  S.   Thus 

^  Y^  7;rT2   CT  (^B2YYY  +  B^Ryyy)  +  ...    and 

Hpy   =   ~    2YYY  ~   2'iYY  "*"••• 
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We  can  see  "by  observation  that  since  the  principle  parts 
of  Hp  and  S  are  proportional  the  last  two  terms  in  the 
first  equation  of  (A1.9)  will  contribute  nothing  to 
the  characteristics  by  virtue  of  the  second  equation. 
The  characteristics  for  the  system  are  thus  found  to  be 

dT^  dY  [dY^  -  ('^L   "^  "^L^^  dT^  ]  =  0  . 

Determining  the  characteristics  of  system  (2.48), 
which  contains  both  magnetic  fields ,  requires  the  same 
type  of  elimination  procedure  used  for  system  (2.52)  only 
now  we  have  three  equations  with  equal  highest  order  terms 
Omitting  the  details  we  simply  state  the  result.   The 
characteristic  equation  for  (2.48)  is 

2     2 

dT^  dY^  [dY^  -  (c+  +  c£  )  dT^]  =  0  . 
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APPENDIX  2.   The  Discrete  Energy  Equation 

We  now  derive  the  energy  equation  (3. 11)  associated 
with  the  difference  equations  used  to  solve  (5.I). 

The  system  (5.1)  written  as  the  difference  equations 
(5.3),  (3.^),  and  (3.10)  is 

62Rn  _  1  Rn(z,+Bn)2  ^  1  Rn^^+^n^   ^0^^ 

2^  j^n     J       ^^n       j^n5 

J  J  J 

r"5V(r"m^)  -  ^  (M^aObOrH^^^  50^^A0R^^ 
J     ^JJ^2^0      J       OJ 

M^  2 

(A2.1)    -  -i  (5°r5)  (A^R^  ))   =   0 

J 

2         >  2 

5"A-(R^  A+B^l)  -  ^  5"(B^A"R^  )   =   0 

J  J  ^  J        J 

^0  n2  2^+1    ^n  2     2 


J  0      J 

J-1    J     J    J- 


-  A+(Z^_^6°R^  +  tP^'^vP      )      =      0 


In  deriving  the  energy  equation  we  follow  the 

procedure  used  to  obtain  (2.54).   Thus,  the  first  equation 

0  n 
in  (A2.1)  is  multiplied  by  5  R.,  the  second  equation  is 

J 

multiplied  by  -  -^  m"?  ,  and  the  third  multiplied  by  -  -^  B . , 
The  term  by  term  results  of  these  multiplications  are 
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=      1    [   I   (5-r5+1)2    -    I    (6-R-)2] 


^  r'^5°R^(A+B^)2   =    _      1      5-(R^+^A+B^+^-    R^A+B^)^ 


^2    "J       "3  '         J 


4  A" 


_      1      A-B^5-[(A+B^+1+A+B^)R^+1r^] 
2  ^^  J  J  J       J         J 


^n,+„nx2       /„n-l^+. 


n-1 


+  ^   [(r5+Vb^+1)2+    (Rn^+B^)2    _    (R;.^A-Bp^-    (r;.^--^A-B  .  1  ^  ] 


1     R^^fiOpn   +  n' 
27^  Rj5    R^.A  Bj. 


^  ^''^^< 


50r^  2    ^      2 

2  AR^  J            <^ 

J  2 

S^R^  i-  = 


1        A  +  Zr^n 


T 


A 


•  O^n 


n^O^n 


J 


J- 


r-7-Li  7  O  Q 


R. 

J 


R. 


2i5-  5°R^ 

J 


R 


+ 


2k 


R^+l-R^-1 


[21 


R^ 


n 


3 


-) 


k:       r.        r 

R^+l    -R^-1 

_^___J 

(R^+1    Rn-1)2 
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■1    „n-,nc- 0.2/^11, 


--^2    r5j55^a'^(r5m5)    ^   ^  5-[(A-(Rn+|M^+l)    _    A-( 


^2    -    LV"    v-j+i-j+iy    -   "    ^Rj+iRj+i))^ 


J  >J  JJ  r,\^  J+J-J+-LJJ' 


2  A" 


•A 


^'^^5.1^5.1)^-^2-  [(--(h-im^:^))^.  (--(H^.i^.i))^ 


+    (A-(r^+1m^+1))2   +    (^-(j,n^))2_    (A-(r^_^^mJ^^))2 


J         J 


-  (^'(R?:K:I))^-  (A-(R^w?))2-  (a-(r^-1m^-1))2j 


j+i  j+i' 


J  J 


J    J 


4  ^^^6°A°R^^ 
2A      J  J 


2  2  2  M^  ? 

^  M^   A^Rj    50^5   =  ^  5-[M^+1^(aOr--^1    +  A^R^    )  ]+  ^{  ^^ [K^Yq    ) 


-^T^F  t^^'^i^^^^r"-^^'^   ) 


.        5V  2    ^      2 


m^m^-1(aV%  a°r"-i^] 


-^  B^A-5  0(R^   A+B^ 
X2      J  ^    J  J 


2^  B^60(B^A-Rn')    =   ^^  5"  (B^-^^B^"  (  R^+l'+R^ )  )  -     ^     S-lB^i^+V^nf 


J         J 


+  A-    [B^+^B^^A-lR^+l    +   R^^    -    bW1(R^%   R^-1^] 
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Notice  that  we  have  substituted  the  fouith  equation  of 
(A2.1)  directly  into  the  sixth  equation  in  the  first  group. 
Each  term  is  multiplied  by  h  and  the  first  group  is  summed 
from  j  =  0  to  J-1  while  the  second  and  third  are  summed 
from  j  =  1  to  J-1.   We  note  the  following  identities'. 

2 
1  ^n  ^n   oO„n      1  T^n^nj-O^ 

-  \   (B^A-60(R^Vb"))^  =  \   (A+b560(R^Vb^))0 
(A2.2)  +  ^  b550(r"j!,A-b5_,).  1,   B-50(Rf  A  B^) 

A  A 

We  now  add  all  terms  and,  remembering  the  definition  (3.11) 
of  K,   we  obtain 
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(A2.5) 


k +  2T  %-l^j   ^   Rj-1 


ij  bS5°(k5_,A+b5.,)    +  -lJ(R5M5.R5_,^^.^6°A-(RX) 


•Ay         "-o 


+  -rr  a 


(z';_.60R^+z':f5°R^  , )  =  \{&-^M)-v  -^(^  +;;^^)  +  ^ 


^J-1      '^J^  -J      '^J-l 


where 


2  2         2 


1 


+1    (5-(B^B^+1),A-r"    )^ 


771    =      \    (M^^A0(R^6V))^-   I    (M"^A05V^^ 


^    (6-(M^+VaO(r"+1%  R^^).1)i+  |(m".6°(M^aV^)^ 


^    =      \    (5-(R^+1a+B^+1-    r'^A+b'^)^1)^-(aV,60(R^    A+B^))-j_ 


+  I    (A+B^,5-((A+b"+1+  A+b")R^+^R^)i 


^=      1    (5-[(A+(R^+V+1)-    A+(R^M^))2],1) 


0 


+  H    (S'LlA-lR'^+V+l)    -   A-(r''m"))2].1)i 


^  = 


2k    ^^^    '^  3 


K 


n' 


2  2 

(r^+1r^-1)2 


and   Y.    =   jh, 

J 


-8> 


^^  '-Vj 


There    remains    only   to   sho'w   that     ^  ,    M,   ^  ,  TPt , 

~        2 

and  ^  are    0(k  )  to  arrive  at  (3.12).   This  Is^  however, 

easily  done,  for  by  simple  manipulation  It  Is  found  that 
these  quantities  can  be  written  as 

S     -    -   ^    {((A+B^)^5-(5-R^+l)2)^-(5-(5-A+B^+l)^R^^^ 

-  ((5-A+B^+l)^5-R^+l^)^+  ((6-A+B^)2,B-R^^,} 
(A2.4)  ^  ^ 

Z?^,  ^   |-{5-(B-A+(R^+1m"+1))2,1)q+  (5-(6-A-(r"+^M^+^))2,1)-^} 

<f  -  -  ^^-  ( ^P-^— ^  .  2R^+1r^-1(B+R^)(5+R^-1) 

j^n-1  ^n-^j^n+l 

[5+R^+  5+R^-1-  R^B^R^]  -  r'^[R^^5-(5+R^)2 


(5+r"-1^(6+R^)2  _  (5+R^^)(5+r'^-1)2]) 


1 


From  (A2.4)  we  see  that  If  R,  B  and  M  possess  sufficient 
bounded  derivatives  then  £>  ,  /K  ,  £>  ,  M  ,    and  (^  will 
all  be  O(k^) . 
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APPENDIX  3.   Boundary  Conditions 

Because  of  the  mixed  hyperbolic-parabolic  nature  of 
the  system  (2.48)  the  number  of  boundary  conditions  one 
would  expect  to  be  necessary  to  determine  uniquely  the 
solution  is  given  by  counting  characteristics  taking  into 
account  infinite  characteristic  speeds.   The  inclusion  of 
pressure  contributes  an  inward  pointing  characteristic  with 


slope  dY/dT  =  -  ycf     +  Cj   .   At  the  inner  boundary  R  =  0 
this  speed  is  seen  to  vanish  so  the  presence  of  pressure 
necessitates  specifying  one  additional  boundary  condition 
at  the  outer  edge. 

From  a  physical  viewpoint  this  condition  should 
reflect  the  fact  that  as  the  plasma  moves  inward  a  large 
pressure  gradient  develops  which  smooths  out  the  edge  and 
prevents  the  formation  of  a  vacuum.   One  way  to  handle  this 
is  to  impose  the  condition  that  the  gas  "sticks"  to  the  wall. 
Numerically,  this  means  that  the  value  of  R  at  the  last 
mesh  point  remains  fixed  at  its  initial  value.   That  is 

(A5.1)  r5  =  Rj  . 

There  are,  however,  some  serious  drawbacks  to  this  condition. 
One  is  that  as  the  density,  and  hence  the  pressure,  becomes 
small  the  problem  becomes  badly  behaved  at  the  boundary  since 
the  existence  of  this  condition  depended  on  the  presence  of 
pressure.   Secondly,  the  large  gradients  in  density  which 
develop  as  the  plasma  is  driven  inward  give  rise  to  numerical 
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instabilities  which  can  only  "be  eliminated  by  taking  a 
prohibitively  small  ratio  of  k/h  in  the  difference  method. 
For  these  reasons  we  do  not  adopt  this  condition.   Instead, 
we  assume  that  the  ion  pressure  effects  at  the  boundary 
can  be  neglected  compared  to  the  force  produced  by  the 
compressing  magnetic  fields.   We  thus  postulate 

(A3. 2)   (RP^Jv  -  ^^2      =       ^—"^h    -   ^  =  0 
±±    X  1    ^^  2YRy^   ^     R-^ 

at  the  boundary.   A  sufficient  condition  for  (A3.2)j  which 
we  use  as  the  boundary  condition  when  a  ,/  0  is 

(A3. 3)  Ny  =  0  . 

The  electron  pressure,  by  virtue  of  its  dependence  on  the 

magnetic  fields  still  contributes  a  force  at  the  boundary. 

However,  as  we  see  from  the  characteristic  equation,  the 

two  pressures  contribute  only  one  additional  characteristic 

so  no  new  boundary  condition  need  be  imposed  when  the 

electron  pressure  is  added.   Since  N  =  2/(R  )„   ,  (A3. 3) 

can  be  written 

2      2 
(A3. 4)  R^   -  2Rj_^  +  Rj_2  =  0  . 

A  heuristic  indication  that  (A3.^)  is  a  reasonable 
boundary  condition  to  Impose,  is  given  by  examining  the 
pressure-free  problem.   In  this  case,  because  of  the  explicit 
presence  of  the  variable  R,  it  is  not  possible  to  solve 
for  Rat  the  last  mesh  point  while  maintaining  the  desired 
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accuracy  of  the  difference  method.   This  arises  because  of 

the  centering  of  the  discrete  variables  and  is  absent  in  a 

plane  geometry.   See  Morton  (I5).   To  determine  R^  in  our 

difference  scheme  we  extrapolate  to  obtain 

2       2      2 
(A3. 5)         ^J   =  2  R^-l  ~  ^J-2  +  °^^^  • 

p 

Better  accuracy  is  obtained  by  extrapolating  on  R  rather 

than  R  because  R  behaves  like  ^/£ .      However,  (A5.5)  is  exactly 
condition  [k^.h)    we  imposed  for  the  problem  in  which  pressure 
was  present.   Thus,  the  condition  (A3.4)  plays  two  roles 
although,  as  a  practical  matter,  its  roles  are  indistinguish- 
able.  One  advantage  to  this  is  that  by  not  introducing  new 
conditions  into  the  problem  we  .can  easily  observe  the  effects 
of  pressure  on  the  disturbances  which  develop  in  the  fluid. 
Of  course,  at  the  same  time  we  lose  all  information  about 
the  boundary  layer  at  the  edge. 

Morton  (I6)  uses  a  similar  boundary  condition  at 
the  plasma-vacuum  interface  in  the  plane  case. 
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